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1.  Intro  due  tion 


Since  spirit  levelling  cannot  be  used  across  the  oceans,  connecting  con¬ 
tinental  vertical  networks  has  long  been  a  challenge  for  both  oceanographers  and 
geodesists.  Among  the  former,  Cartwright  (1963)  calculated  a  tie  between  the 
British  and  the  European  nets  across  the  English  channel.  Because  the  oceanog¬ 
raphic  method  requires  a  knowledge  of  currents  that  is  not  available  for  larger 
bodies  of  water,  the  geodesist  Erik  Tengstrom  (1965)  tried  using  gravimetry  and 
deflections  of  the  vertical  to  compute  a  geoidal  profile  through  the  Eastern  Medi¬ 
terranean,  from  Athens  to  Alexandria.  His  results  suffered  from  lack  of  data. 
Lelgemann  (1976)  proposed  unifying  vertical  datums  by  means  of  gravimetry, 
levelling,  and  very  accurate  position  determinations,  as  those  expected  by  the 
proponents  of  lunar  laser  ranging  (Sllverber  et  al. ,  1976).  The  late  R.  S. 

Mather  (1978)  considered  the  possibilities  open  for  datum  unification  by  constant 
improvements  in  gravity  field  models  (notably  Goddard's  GEM  series)  and  the 
large  amounts  of  information  provided  by  the  altimeter  satellites  whose  proto¬ 
type  has  been  GEOS-3. 

Dynamic  effects  create  a  "sea  surface  topography",  or  departure  of  the 
mean  sea  surface  from  a  true  equipotential.  This  topography  is  not  very  well 
known  at  present,  and  has  an  estimated  r.  m.  s.  value  of  about  1  m.  Tide 
gauges  determine  the  position  of  mean  sea  level  at  their  locations,  so  the  error 
made  by  assuming  that  their  mean  sea  level  marks  are  on  the  same  equipotential 
surface  (geoid)  is  of  the  order  of  /2  x  (r.  m.s.  of  the  stationary  sea  surface 
topography)  1.5kgalm.  Without  any  additional  work,  Nature  provides  a  world 
"levelling  net"  of  1.  5  kgal  m  accuracy. 

Would  it  be  possible  to  obtain  better  transoceanic  links  using  the  various 
forms  of  geodetic  data  that  are  available  at  present,  or  are  likely  to  become 
available  in  the  near  future?  Could  it  be  feasible,  with  such  data,  to  establish 
benchmarks  for  levelling  inside  continents,  rather  like  inland  "tide  gauges", 
whose  potential  differences  are  known  so  well  that  they  can  be  used  to  constrain 
the  adjustment  of  the  net  to  reduce  distortion?  If  so,  in  a  more  distant  future, 
similar  benchmarks  could  be  used  to  survey  other  components  of  the  Solar  System, 
where  only  the  Earth  has  any  significant  amount  of  free  surface  water. 

In  this  work  the  reader  will  not  find  more  than  a  passing  reference  to  the  geoid, 
a  notion  that  appears  inseparable  from  that  of  vertical  height  and  of  vertical  datum. 
Since  the  geoid  has  been  regarded  as  the  natural  universal  datum  by  geodesists,  a 
few  words  of  explanation  are  due.  The  reason  for  its  omission  here  is  that,  however 
useful  otherwise,  the  geoid  is  not  essential  to  the  setting  up  of  a  vertical  network, 
at  least  in  theory.  Such  network,  ultimately,  is  a  set  of  potential  differences  es¬ 
timated  among  the  points  that  form  the  net,  in  particular  the  primary  points  or 
benchmarks.  These  potential  differences  do  not  convey  information  on  the  abso¬ 
lute  potential  of  the  gravity  field,  so  they  can  be  referred  to  any  number  of  level 
surfaces,  and  not  exclusively  to  one.  For  the  same  reason  their  meaning  is  not 


dependent  on  which  surface  is  selected,  and  it  remains  intact  even  if  no  surface 
is  selected.  The  reason  why  the  geoid  is  so  ubiquitous  in  the  literature  on 
levelling  may  be  the  complete  reliance  that  levelling  has  had  on  tide  gauges, 
idealized  as  points  on  the  geoid.  As  long  as  gauges  play  a  basic  role,  the 
concept  of  reference  surface  is  relevant  to  levelling. 

Here,  instead  of  heights  above  an  equipotential  surface,  we  are  going  to 
consider  distances  to  the  center  of  mass  of  the  Earth,  or  to  a  reference  ellipsoidal 
centered  on  this  point.  Such  approach  is  not  unreasonable  today,  when  new  po- 
sitioning  techniques  are  being  developed  that  promise  accuracies  far  better  than 
those  available  in  the  past.  Methods  based  on  the  Global  Positioning  System 
satellites,  on  Lageos,  and  on  portable  interferometric  and  lunar  laser  ranging 
stations,  are  expected  to  achieve  near  decimeter  accuracy  in  relative  position, 
over  continental  distances. 

In  recent  years,  the  use  of  artificial  satellites  has  changed  many  aspects 
of  geodesy.  Space  techniques  for  obtaining  position  fixes  and  models  of  the 
gravity  field  are  in  constant  development,  and  both  the  quality  and  the  quantity 
of  the  data  provided  by  spacecraft  are  increasing.  This  work  shall  explore 
how  these  advances  may  affect  levelling.  Next  paragraph,  to  begin  with,  intro¬ 
duces  a  basic  idea*  a  World  Vertical  Network  established  without  recourse  to 
any  reference  surface,  or  geoid.  In  a  way,  such  network  is  the  datum. 

1.1.  Definition  of  World  Vertical  Network 

The  World  Vertical  Network  (WVN)  is  a  set  of  estimated  potential  differences 
among  benchmarks  situated  in  various  continents. 

A  network  of  potential  differences  can  be  translated  immediately  into  a  variety 
of  height  systems  such  as  those  described  by  Krakiwsky  and  Mueller  (1965).  Po¬ 
tential  differences  are  intrinsic  to  the  set  of  benchmarks  selected,  and  are  indepen¬ 
dent  on  the  choice  of  "geoid",  or  on  the  precise  knowledge  of  the  zero  harmonic 
of  the  Earth's  potential,  present  estimates  of  which  have  an  uncertainty  of  some 
3  kgal  m.  Existing  regional  networks  can  be  tied  to  the  closest  benchmarks  to 
create  a  dense,  unified  global  levelling  net.  Any  benchmark  potential  can  be  used 
to  reference  all  points  tied  to  the  net.  If  so  desired,  the  level  surface  through 
this  arbitrary  point  may  be  regarded  as  a  "geoid". 


2.  Method  of  Approach 

If  we  had  a  perfect  model  of  the  gravity  field  and  exact  position  fixes  in 
geocentric  coordinates  at  two  points  on  the  Earth's  surface,  then  we  could  use 
this  Information  to  find  the  potential  of  each  point  and,  from  this,  their  potential 
difference.  Repeating  this  process  for  all  possible  pairs  of  points  out  of  a  given 
set  of  benchmarks,  the  end  result  would  be  an  exact  WVN.  Unfortunately,  models 


and  fixes  are  never  perfect,  so  the  potential  differences  must  have  some  errors. 
To  reduce  these  errors,  we  could  combine  the  field  model,  which  being  finite 
cannot  contain  information  above  certain  spatial  frequencies,  with  additional 
data  such  as  gravimetry,  rich  in  high  frequencies,  particularly  from  the  vicinity 
of  the  benchmarks. 


2.1.  Formulation  of  the  Problem 

If  V  is  the  gravitational  potential  due  to  the  mass  of  the  Earth  and  external 
to  its  surface,  and  U  is  a  reference  potential  defined  by  a  spherical  harmonic's 
model* 


U (<p,\  , r)  =  - [1  +  £  (a/r)BPal(sin<fi)  fCn,cosmA+Sn,sinmA}]  (2*  *) 

r  a  =a 

where:  Pn.  hilly  normalized  Legendre  function  of  the  first  kind,  degree  n  and 

order  m  ; 

r,<fl,A  geocentric  distance,  latitude,  and  longitude; 

G  universal  gravitational  constant; 

M  mass  of  the  Earth; 
a  mean  equatorial  radius  of  the  Earth; 

Ca„Sn.  normalized  spherical  harmonic  coefficients; 

N  maximum  degree  and  order  for  terms  present  in  the  model; 
then  the  disturbing  potential  T  at  a  point  P  of  geocentric  coordinates  r(.,<pP,Ai>  is 


T  (  P )  =  V(P)  -  U  (P)  (2.2) 


The  gravity  potential  of  the  Earth  is 


W(P)  »  V(P)  +  <P  (P)  (2.3) 

where  <o(P)  =  £  oc3  r,3 cos3 corresponds  to  the  rotational  potential,  u:  being 
the  angular  velocity  of  the  planet  about  its  spin  axis.  The  potential  difference  be¬ 
tween  two  points  such  as  P  and  Q  is,  therefore, 


iW(P,Q)  =  U(P)  +  T(P)  +  <P(P)  -  U(Q)  -  T (Q)  -  <P(Q)  (2.4) 

With  both  P  and  Q  on  the  Earth's  surface,  the  uncertainties  in  the  calculated 
values  of  <P(P)  or  <fi(Q) ,  due  to  errors  in  the  known  positions  of  P  and  Q  ,  will  be 
thousands  of  times  smaller  than  those  arising  in  the  determination  of  U  (P),  U(Q), 
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T(P)  and  T(Q) .  Consequently,  only  errors  in  the  computed  values  in  the  right 
hand  side  of 


iW(P,Q)  -IflP)  -*(Q)]  =  V(P)  -  V (Q) 

-U(P)  +  T(P)  -  (U(Q)  +  T(Q)]  (2.5) 

shall  be  included  in  the  error  analysis  that  constitutes  the  major  part  of  what 
follows. 


2.2.  Estimating  T  by  Least  Squares  Collocation 

The  use  of  linear  regression  for  predicting  and  filtering  geodetic  data 
appears  to  have  been  first  proposed  by  Kaula  (1959).  Further  developed  by 
Moritz  and  others,  this  approach  has  become  a  familiar  technique  that  has 
shown  its  value  in  many  applications.  Least  squares  collocation,  as  geodesists 
call  it,  provides  a  way  of  combining  all  relevant  data  into  estimates  of  unobserved 
variables  (minimum  variance  prediction),  or  into  more  reliable  estimates  of 
those  actually  observed  (minimum  variance  filtering).  For  further  information  on 
this  method,  see  Moritz  (1972). 

If  T  is  to  be  estimated  at  point  P,  then  the  linear,  unbiased,  minimum 
variance  estimator  of  T  (P)  is 


T (P)  -  f1 '  d  =  fT(z  +n)  (2.6) 

where 

l  =  (C  n  +  D)-1  C  rz  (2.7) 


is  the  optimal  estimator  vector,  and 


A 

T 


d  =  z_  +  n 

is 

the 

z 

is 

the 

n 

is 

the 

Crz 

-  MfTzr] 

is 

the 

Cu 

=  M{££T  } 

is 

the 

D 

=  M  [n_nT } 

is 

the 

is  the  estimated  disturbing  potential,  a  scalar; 

Nj  vector  of  measurements,  or  data  vector; 

vector  of  signal  component  in  the  measurements; 

Ni  vector  of  the  noise  component  in  the  measurements 
lxNi  covariance  matrix  (a  row  vector)  of  T  and  z; 
Nj  xNj  covariance  matrix  of  the  signal  z  ; 

Mi  xN4  covariance  matrix  of  the  noise  n. 


The  operator  M  {  1  represents  some  kind  of  average.  D  is  supposed  to  be 
diagonal,  because  the  noise  is  not  correlated  from  measurement  to  measurement 
(white  noise).  Furthermore,  these  assumptions  apply: 


.UtiiJIPPj  jpii 


and 


M  [_z  }  =  M  {n}  =  M  [d]  =  £  (a  null  vector) 

M  [z  nr}  ,  M  f  T  nT  ]  are  both  null  matrices. 


More  generally,  we  could  be  asked  to  obtain  N,  estimates  £  (N,  vector  of 
estimates)  from  d ,  using  an  estimator  matrix  F  such  that 


A 

s 


Ft  d 


(2.6)* 


minimizes  the  mean  square  values  of  the  components  of  the  error  vector 


e 


s  - 


A 

s 


(s  is  the  Ns  vector  of  true  values  of  s).  The  variance-covariance  matrix  of 
these  errors  is 


E  =  M(eer}  =  Mf(a  -  Frd)(s  -  FTd)T}  =  C„- FrC,Tt-C.tF  +  Fr(C„  +  D)  F  (2.8) 

where  C„  =  M{s,sT }  is  a  N,x  N,  matrix,  and  C,z  =  M  [s  1}  is  a  N.x  Nt 
matrix.  Since  the  elements  in  the  main  diagonal  of  E  are  either  positive  or 
zero,  minimizing  each  one  of  the  mean  square  errors  is  the  same  as  minimizing 
their  sum,  the  trace  of  E  (tr(E)).  Accordingly  (see  for  instance,  (Rao,  1973)), 

§-§JiEl£)  =  cjj  +  (c„  +  D)  F  =  fif  (a  null  matrix) 


or 


(Czz+D)  F  =  C,7. 

(2.9) 

and,  finally, 

F  =  (Cu  +  D)"1  C.\ 

(2.10) 

Replacing  (2,10)  in  (2,8)  we  get 

E  =  C„  -  C a 2  (Cjj  +  D)  Cft 

(2. 11) 

In  the  special  case  where  s  is  the  scalar  T  ,  we  get  (2. 7).  The  equations  in  the 
system  (2.9)  are  known  as  "normal  equations";  some  people  prefer  to  call  them 
"Wiener-Hopf  equations"  because  they  bear  a  formal  resemblance  to  the  basic 
integral  equation  of  linear,  invariant,  minimum  variance  filtering  in  the  time  domain. 
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In  geodetic  applications,  M  f  }  is  an  average  on  rotations  of  some  sort. 

If  all  possible  rotations  about  the  origin  (center  of  mass)  are  included,  then  the 
covariance  function  c^  of  a  function  h  of  r ,<p  and  A 

M{h(P)  h(Q)}  *  chh(P,Q)  (2.12) 

depends  only  on  the  spherical  distance 

0*,  =  cos-1  [sin  <A>  sin  <4;  +  cos  cop  cos  <qj  cos  (A*  -  Aq )] 


and  on  the  geocentric  distances  r?  and  rQ .  If  both  P  and  Q  are  on  the  same 
sphere  ( r,,  =  rp),  then  chS  depends  on  (Dpq  alone.  For  this  reason  this  type  of 
covariance  is  known  as  isotropic,  and  the  operator  M{  }  is  then  called  the  iso¬ 
tropic  average  opera  tor7]  The  elements  of  F  depend  on  those  of  matrices  C.i 
and  C«,  and  these  elements  are,  in  turn,  values  of  the  covariance  functions 

c»i(P,Q)  =  Mfs  (P)  z  (Q)}  and  c„(p,Q)  =  m[z(P)  z(Q)}. 

A  choice  of  M[  }  determines  those  functions,  their  values,  and,  ultimately,  the 
optimal  estimator  matrix  F .  Rummel  and  Schwarz  (1977)  have  discussed  different 
types  of  averages  and  covariance  functions.  From  these  considerations  it  is  clear 
that  the  optimal  estimator  is  not  unique,  but  it  depends  on  what  average  we  choose. 
The  "easiest"  choice  is  the  isotropic  average,  because  of  the  simplicity  of  the 
corresponding  covariance  function.  This  function  can  be  expanded  as  a  series  of 
Legendre  polynomials 

00 

Cuu(P.Q)  =  y  (2n  +  l)(_lLN!  6  ,  Pa(cosv*))  (2.13) 

n"sto  r,/ 

where  6uu,n  is  the  nth  degree  variance  of  the  spherical  harmonic  coefficients  of 
u( <P,A , r ) : 

n 

„  =  Y  (  Cu3,n.  +  C- )  ( 2n  +  1  )_1  (2-  14> 

•  =  0 

The  covariance  between  two  functions  u  and  v  is 

1  There  is  a  small  problem  with  Mfn  nT},  because  the  measurements'  "noise"  is 
notan  ordinary  function  of  <0,  A  and  r,  but  a  stochastic  process.  However,  it 
can  be  manipulated  as  if  it  were  such  a  function.  For  this,  see  the  discussion  by 
Balmino  (1978). 


(2. 15) 


/  a3  Y 

<V,„(P,Q>  -  t  ( 2n  +  1)  VrTrTV  &^.n  PB(cos  <&*,) 

0*0  * 

where  * 

a 

5iiv,n  =  ^  (Cu,BSC»,ni  +  Su,„  t  S»  ,Bi )  (2  n  +  1 )  (2.16) 

laO 

To  understand  in  what  sense  the  estimator  is  "optimal",  imagine  some 
pattern  of  measurement  points  and  estimation  points.  The  whole  pattern  is 
subject,  in  succession,  to  ail  possible  rotations.1  Before  each  rotation, 
measurements  are  taken  and  all  estimates  are  made  at  their  respective  points, 
and  the  squares  of  the  estimation  errors  are  found,  somehow.  This  is  repeated 
over  and  over  again,  and  running  averages  of  the  errors  squared  are  kept.  In 
the  limit,  these  averages  will  tend  to  values  that  satisfy  (2.8);  if  F  is  optimal, 
they  will  also  satisfy  (2. 11)  and  will  be  smaller  (or  not  worse)  than  for  any 
other  choice  of  F.  Also,  in  the  limit,  we  would  have  covered  the  whole  Earth 
with  estimates,  which  is  why  such  mean  squared  values  and  their  square  roots 
(r.  m.  s.  values)  are  called  global.  In  practice  we  are  always  concerned  with  a 
finite,  even  a  small  number  of  estimates  at  isolated  locations,  and  we  are 
interested  in  the  actual  errors  of  those  estimates,  not  "some  global  measure". 

The  practical  meaning  of  the  latter  is,  therefore,  a  matter  of  interpretation.  If 
signal  and  noise  have  near  Gaussian  distributions,  then  the  errors  (which,  accor¬ 
ding  to  (2.6),  (2.6)*,  are  linear  transformations  of  both)  will  also  be  near 
Gaussian.  In  such  cases  the  global  values  are  related  to  the  actual  errors  by  the 
usual  "one  sigma"  and  "three  sigma"  rules,  giving  an  indication  of  their  likely 
sizes.  Rapp  (1978a)  has  shown  that  a  world-wide  data  set  of  38406  l°x  1°  mean 
gravity  anomalies,  compiled  at  The  Ohio  State  University,  has  a  nearly  Gaussian 
distribution.  Gravity  anomalies  are  the  main  type  of  data  considered  in  this 
report  for  predicting  T. 

The  probability  distribution  of  the  data  does  not  characterize  it  enough, 
however,  because  all  the  large  values  could  be  concentrated  in  a  few  "rough" 
areas,  the  rest  of  the  world  being  "smooth",  with  smaller  values.  The  errors 
are  likely  to  repeat  this  pattern  (see  Appendix  A)  so,  if  estimates  are  made  in  a 
"rough"  region,  the  global  r.  m.  s.  may  give  an  over-optimistic  indication  of  the 
actual  size  of  the  errors.  This  quality  of  the  data  being  "evenly  behaved",  so  that 
there  are  no  zones  that  are  highly  idiosyncratic,  is  known  as  stationarity.  It  is  a 
rather  elusive  quality,  but  very  important  to  the  use  of  global,  isotropic  covariances. 
How  stationary  is  the  Earth?  We  know  that  trenches  and  ridges  in  the  ocean  floor 
produce  strong  localized  features  in  the  gravity  field,  set  off  against  comparatively 
smooth  surrounding  areas.  Mountainous  regions  in  land  also  tend  to  be  "rough"; 
however,  there  are  very  flat  regions,  such  as  the  Nullarbor  plain  in  S.  W.  Australia, 
where  the  field  presents  strong  local  anomalies. 


1  Not  only  rotations,  but  more  generally  all  orthogonal  transformations  can  be  in¬ 
cluded  (i.  e. ,  rotations,  reflections  and  various  symmetries),  the  result  being 
precisely  the  same  average  values  as  with  rotations  alone. 


Regardless  of  the  significance  of  the  results,  getting  them  can  present  dif¬ 
ficulties.  We  have  to  form  and  invert  a  matrix  whos'»  dimension  is  that  of  the 
data  vector,  so,  if  many  measurements  are  involved,  this  two  operations  become 
quite  large.  Not  only  the  computer  time  involved,  but  also  the  accumulation  of 
rounding  errors  can  escalate  dramatically.  Paradoxically,  the  more  data  are 
used,  the  better  the  results  (in  theory),  but  also  the  harder  to  get  and  the  more 
unreliable.  This  is  further  complicated  by  the  fact  that,  for  close  spacings  of 
data,  the  normal  equations  can  become  very  ill-conditioned.  A  special  technique, 
presented  in  Section  4,  has  been  developed  by  the  author  to  overcome  these 
problems  in  the  case  at  hand. 

On  the  positive  side  we  must  consider:  the  possibility  of  using  mixed  data 
sets,  so  d  may  consist  of  gravity  measurements,  satellite  altimetry,  deflections 
of  the  vertical,  and  even  levelling;  the  ability  to  provide  more  than  one  optimal 
estimate  at  the  same  time;  the  simplicity  and  elegance  of  the  theory.  Another 
good  aspect  of  collocation  is  that  the  covariance  functions  needed  to  set  up  C,z 
and  ( C«  +  D)  do  not  have  to  be  known  with  great  accuracy.  This  is  born  out 
by  the  results  presented  in  Section  5,  where  the  same  problem  has  been  solved 
using  somewhat  different  covariances.  This  is  fortunate,  as  we  can  never  gather 
sufficient  data  to  obtain  an  exact  empirical  covariance,  because  to  know  such 
function  is  equivalent  to  knowing  the  whole  field  exactly  ( thus  making  estimation 
unnecessary). 


X 


Figure  2. 1.  The  circles  represent  spherical  caps  within  which  gravity 
anomalies  with  respect  to  a  reference  model  have  been 
measured.  The  wandering  lines  are  levelling  traverses. 

T  is  estimated  at  the  center  of  each  cap.  The  differences 
in  the  geocentric  distances  of  these  center  points  are 
known  to  decimeter  accuracy.  BMA  and  BMB  are  two 
benchmarks  of  the  WVN. 


Figure  2. 1  shows  the  basic  data  arrangement  to  be  studied.  While  for  the 
simulations  of  Section  5  all  caps  are  supposed  to  have  the  same  size,  to  reduce 
computing,  this  is  by  no  means  essential.  Other  kinds  of  information  (such  as 
satellite  altimetry)  could  be  included  among  the  data  (see,  for  instance,  paragraph 
5.3),  though  only  the  types  indicated  in  Figure  2.1  shall  be  considered  here. 
Anomalies  are  determined  at  the  Earth's  surface,  somewhat  in  the  manner  of 
Molodenskii.  Details  are  given  in  Section  3. 


2. 4.  Adjustment  Theory 

Consider  a  cap  center  Pt  in  zone  "A"  (Fig.  2. 1)  and  another  Qj  in  zone 
"B".  The  potential  difference  between  benchmarks  BMA  and  BMB  is 

AW(BMA.BMB)  =  U(Pi)+T(P$)  +  AW^ (BMA,  P,)  +  <P(Pi)  -  U(Qj) 

-  T(Qj)  -  aW£(BMB,QJ)-0(Q4)-V1J  (2.17) 

where 

V tj  =  *AW^(BMA,P,)  -  fAWj (BMB.Qj)  +  eU(Pi)  -  CU(Qj) 

+  eT(Pi)  -  eT(Qj) 


is  the  sum  of  the  error  in  levelling  e  AW^  ,  the  error  in  the  potential  according  to 
the  reference  model  e  U  ,  ar.d  the  error  in  the  estimated  disturbing  potential  cT.1 
(U  is  influenced  both  by  errors  in  the  model  and  errors  in  the  coordinates  of  the 
p, ,  Q. .  Expression  (2. 17)  can  be  regarded  as  an  observation  equation  with  the 
potential  difference  AW  as  the  only  unknown.  Each  pair  of  caps  (Pi,Qj)  pro¬ 
vides  an  equation  of  this  kind,  so  a  redundant  system  can  be  set  up 

a  AW(BMA.BMB)  =  £  +  v  (2.18) 

where  £  is  the  vector  of  "observed"  potential  differences,  v  is  the  vector  or 
residuals,  and  a  is  a  vector  with  all  components  equal  to  1 ;  the  design  "matrix" 
of  this  particular  system.  All  three  vectors  have  for  dimension  the  number 
of  equations.  This  number  is  restricted  by  the  following  considerations:  the 
centers  of  the  caps,  paired  in  the  same  way  as  the  caps,  should  not  form  a  closed 
loop,  as  shown  in  Figure  2.2  by  a  broken  line.  Otherwise,  because  the  "observed" 
values  for  the  pairs  in  the  loop  always  add  up  to  zero,  i.  e. ,  are  linearly  dependent, 
the  a  priori  variance-covariance  matrix  of  the  "observed"  values  must  be  singular. 

1  As  explained  in  Section  2.1,  errors  in  <fl(Pi)  and  <P(Q.)  due  to  errors  in  the  co¬ 
ordinates  of  Pi ,  Qj ,  and  in  the  rotation  rate  on,  are  considered  to  be  negligible  here. 
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The  inverse  of  this  matrix  is  needed  for  the  adjustment  of  AW,  so  this  limits  the 
choice  of  pairs  of  caps  to  those  not  forming  loops.  The  number  of  such  pairs  is 


Qi  Qa  Q3  Q*  "B" 


Figure  2.2.  The  solid  and  broken  lines  identify  caps 
that  have  been  "paired".  The  broken  line 
shows  a  selection  containing  a  loop  (not 
permitted).  The  maximum  number  of 
equations  (permitted)  =  number  of  caps-1. 


The  accuracy  of  AW(BMA.BMB)  after  the  adjustment  can  be  computed  from  the 
following  formula 


*  (iT  V“1a}’i  =  (E  Z  (VU*  (N.  *  no.  of  equations  (2.19) 

'  *  «  no.  of  caps  -  1  ) 

where  V,  the  variance -covariance  matrix  of  the  data,  is 


V  "  V€Aw£+  VCA.  +  v«At  (2.20) 

In  this  expression, 

V( is  the  variance-covariance  matrix  of  the  levelling  errors; 

V€£g  is  the  variance -covariance  matrix  of  the  errors  in  U(Pi)  -U(Qj) ; 
is  the  variance -covariance  matrix  of  the  errors  in  ’t(Pi)  -^(Qj) . 


The  uncertainties  in  AW^  and  AU  depend  on  those  of  the  data  that  give  them 
origin;  eA^1  depends  also  on  the  way  T  is  estimated  from  the  gravity  anomalies. 


IT 


3.  Characteristics  of  the  Data 


In  this  section  we  shall  consider  those  characteristics  of  the  data  that  affect 
the  accuracy  of  potential  differences  adjusted  according  to  the  method  explained  in 
the  previous  paragraph.  The  various  types  of  data  involved  are:  position  fixes, 
the  coefficients  of  a  reference  gravity  field  model,  gravity  anomalies,  and  levelling 
traverses. 


3.1.  Position  Fixes 

To  compute  the  reference  potential  U  at  the  center  of  every  cap,  it  is 
necessary  to  know  the  coordinates  of  these  points.  The  reference  model  contains 
terms  below  degree  20  or  30  (in  this  study)  so  the  smallest  detail  it  can  show  is 
of  the  order  of  1000  km.  U  is,  therefore,  a  smooth  function  of  latitude  and  longi¬ 
tude,  and  cannot  be  substantially  affected  by  horizontal  position  errors  of  the  order 
of  less  than  2  m,  which  is  the  accuracy  that  can  be  obtained  at  present  with  satellite 
Doppler  techniques.  The  reference  potential,  on  the  other  hand,  is  quite  sensitive 
to  vertical  (geocentric  distance)  errors,  at  the  ratio  of  about  1  kgal  m  per 
meter  of  error.  As  we  are  concerned  with  finding  potential  differences, 
the  most  important  errors  are  those  in  relative  vertical  position.  To  be 
precise,  the  vertical  position  of  interest  is  the  distance  to  the  geocentre.  There 
is  little  difference,  however,  between  relative  errors  in  ellipsoidal  heights  and 
relative  errors  in  geocentric  distance,  and  both  can  be  regarded  as  equivalent  here. 
The  absolute  vertical  error  might  contain  a  nearly  constant  bias,  due  to  the  incor¬ 
rect  dimensions  of  the  reference  ellipsoid  and  to  other  systematic  causes  related 
to  the  positioning  method.  This  error  may  be  of  several  meters  without  any  notice¬ 
able  effect  on  the  estimated  potential  differences,  because  it  will  nearly  cancel-out 
when  such  differences  are  between  points  at  the  Earth's  surface.  Therefore,  an 
error  in  the  ellipsoid  of  the  order  of  2  m,  which  is  the  present  level  of  accuracy, 
can  be  disregarded. 

While  the  relative  vertical  position  error  is  the  one  that  matters,  we  stlU 
have  to  know  the  absolute  geocentric  distance  to  compute  U .  This  can  be  done, 

'  essentially,  in  two  ways  j 
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a)  find  the  absolute  position  of  each  point  separately; 


b)  find  the  absolute  position  of  one  point,  and  then  obtain  the  relative 
position  ef  the  other  points  with  respect  to  this  one. 

In  each  case,  the  relative  position  errors  will  vary,  the  choice  being  always 
the  alternative  that  gives  the  smallest  errors. 

This  study  is  more  concerned  with  forthcoming  developments  than  with  the 
present  state  of  affairs:  Anderle  (1978)  has  estimated  that  the  Global  Positioning 
System  currently  being  deployed  could  provide,  when  all  the  satellites  are  oper¬ 
ational  in  the  mid-Eightles,  relative  positions  with  errors  of  less  than  0.1  m. 
This  accuracy  should  be  possible  between  stations  thousands  of  kilometers  apart, 
in  all  three  coordinates,  after  less  than  one  day  of  constant  observation  of  the 
satellites.  Estimates  for  absolute  position  determinations  from  lunar  ranging 
stations  made  by  Silverberg  et  al.  (1977),  using  mobile  stations  supported  by  a 
network  of  a  few  fixed  ones,  are  also  in  the  decimeter  range.  La  addition  to 
these  two,  a  variety  of  new  positioning  techniques  based  on  satellites  in  high 
orbits  that  carry  laser  reflectors,  mobile  radiointerferometry,  etc. ,  being 
investigated  at  present,  might  provide  even  better  accuracies  in  the  coming 
decade.  Present  measurements  from  Doppler  satellites  have  errors  that  are 
one  order  of  magnitude  worse.  However,  considering  the  progress  made  in  this 
field  over  the  past  decade,  and  the  new  highly  precise  methods  in  the  offing,  it 
is  probably  not  too  optimistic  to  assume  in  this  study  relative  accuracies  as  good 
as  one  decimeter  in  vertical  position. 


3.2.  Reference  Model  of  the  Gravity  Field 


The  coefficients  CB1,SB„  and  the  constant  GM  in  (2. 1)  are  not  exact  values, 
so  the  model  does  not  represent  to  perfection  the  first  N  harmonic  degrees  of 
U  or  T.  The  effect  of  an  incorrect  GM  is,  at  the  present  level  of  accuracy, 
equivalent  to  a  bias  of  about  ±  3  m  in  geocentric  distance  (Lerch  et  al. ,  1978). 

This  error,  being  virtually  constant,  has  a  negligible  effect  on  potential  differences. 
The  existence  of  coefficient  errors 

=  ^a»(tru»)  “ 

»  3n.(«ru.)  -  S,.  <3' 

has  to  be  considered  when  defining  the  disturbing  potential1 


1  The  0  degree  error  cGM/r,  being  almost  constant  on  the  Earth's  surface  has 
no  relevance  to  this  work,  and  has  been  excluded  from  these  formulas. 
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H  ft 


T (P)  =  V(P)  -U(P)  =  (~jy  Pai(sin^)[€CatcosmA,  +  cSn.sinmA,] 


i  =  2  a=0 

38  ft 


7  (“f  Jpn.(slnifi')  [  Ca^cosmAf  +  S„sinmA/>)}  (3.2) 


and  the  gravity  anomaly 


^VT(P)  *y(P)  T(P)  GM  rr/a  Y  V  -  -  - 

Ag(P)  =~~Jr~^+  yfPl  p«<sin^>[eC«cos^  +  fS»*sinmA',] 

'  \  t  ^  a  —  2  isO 

oo  ft 

+  7  f—  J(n-1)Y  Pn,<sin<A>)  [  C„cosmAp  +  3a.sinmAp]j> 

,=°  (3.3) 

The  corresponding  isotropic  covariances,  according  to  expression  (2, 13)  are 

c  rr  (P,  Q)  =ii-23-{)  (2n  +  l)6€nPR(cos^Q)(-— '  +  >  (2n  +  l)6rPa(cos0Pq)(— )  \  (3.4) 

rf  Vr'  ’Vnitf  +  i  ^  V  J 

2  g  N  2  f*  ^  2  w 

cA*&*(P.Q)  »  “^{^(2n+lXn-lf«€aPn(cos^)^)+^(2n+lXn-lf^Pn(cos«w)(jM]-  (3.  5) 

^  n-2  ^  ft  =N  +1  ** 

2  _2  ^  g  co  g 

Cr  A.  (P.  Q)  *  ““"3  2(2n+1)<n'1)6CnP,'(COS,i)pQ)(^)  +  Z(2n+1)(n"1)6:lPa(COS0P9)(^')j  (3‘  6) 


where 


and  6n  s  fiTr,n 


u 

H  *  7  (*CB/  +  €Snsa)/(2n+l) 


In  practice,  neither  6en  nor  6n  (n>  N)  are  known  from  direct  measurements. 
The  6(n  can  be  approximated  by  the  corresponding  diagonal  elements  of  the  a 
posteriori  variance-covariance  matrix  of  the  adjustment  that  produced  the  model. 
The  6n  are  known  to  follow  a  more  or  less  asymptotic  decline  law  such  as 


known  as  "Kaula's  rule",  or  such  as 

,  /GMVV  ax Ri3  (tS _ SaR l _ (&)* )  (3  8) 

°»  "  \  a  A(2n+l)(n-l)(n  +  A)\?y  (2n  +  l)(n-l)(n-2)(n  +  B)Va  /  J  '*  ; 

This  last  expression,  a  "two  terms  law",  has  been  investigated  by  jekeli  (1978), 
using  the  available  information  on  the  power  spectrum  of  gravity  anomalies,  geoid 
undulations,  and  the  horizontal  gradient  of  gravity,  to  find  the  parameters  A,  B, 
ttu  «!)  Bn  Rj  that  give  the  best  fit  to  such  information.  This  type  of  formulas 
has  the  advantage  that  the  covariances  (3.4),  (3.5),  and  (3.6)  can  be  computed  using 
finite  recursions. 


3.3.  Propagation  of  Position  Errors  through  the  Computed  U 


Because  of  the  low  degree  of  the  terms  in  the  expansion  of  U  ,  horizontal 
errors  in  tfl  and  X  are  of  little  importance.  The  vertical  errors  er ,  on  the  other 
hand,  have  a  significant  effect.  If  they  are  small,  we  can  write 


€  AU(P,Q)  =<U(P).CU(Q)-^[|?UB(P)e,p-|?Un(Q)«f,'’  (3.9) 

n  =  0 


where  U„ 

and 


is  the  nth  harmonic  of  U  .  Calling 

AUa(P,Q)  =  U,(P)  -U„(Q) 


we  have 


jeAUj  *  2<n  +  1}  max 
aa(r,0) 


(3. 10) 


where  r  =  min  ( v? ,  rq)  and  ^  a  (  r,0)  is  the  surface  of  the  sphere  o(  r,0 ).  Since 

~  f  i®.D 

for  n  >  0,  it  follows  that 

|eAU0!  >> 


for  n  >  0  .  Calling 

?  ■  i  ( r,  +  er>  +  n,  +  erq) 

we  get 


Uau|  -  leAUo!  “  yr  |crpQ|  ="  Y  |c, J  (3.11) 

where  y  is  the  mean  value  of  gravity  acceleration  on  the  Earth's  surface 
(y  ** 0.9798  Kgal).  The  standard  deviation  of  |eAu|  is 


0  eAu  y  °  €r 

rpq 


3 

l 
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(3. 12) 


If  the  coordinates  of  points  P  and  Q  were  determined  separately,  so  c  r,  and 
c  hj  could  be  considered  uncorrelated,  and  if  or  r*  =  ?crg  =  a  cr  then 


<7  f  4U  *■  y  o  c. 


(3. 13) 


On  the  other  hand,  if  the  difference  in  vertical  position  is  determined  simultaneously 
for  P  and  Q,  then  (3. 12)  applies.  In  any  case,  if  all  geocentric  distances  are 
computed  with  uncorrelated  errors,  except  for  some  constant  bias  that  does  not 
affect  the  potential  differences,  then  matrix  Vf  ^  in  (2.20)  must  be  diagonal, 
each  non-zero  term  being 


vlt  =  c^fhui  (3.14) 

where  Aut  is  the  potential  difference  between  the  centers  of  the  ith  pair  of  caps. 


3.4.  Gravity  Anomalies 


The  theory  used  in  this  work  assumes  that  the  exterior  potential  of  gravi¬ 
tation  is  harmonic.  This  is  not  strictly  correct  at  the  Earth's  surface, 
because  of  the  atmosphere  above  it.  To  avoid  systematic  errors  this 
effect  should  be  discounted  from  the  measured  gravity  values,  and  put  back  on 
the  estimated  potential.  These  atmospheric  corrections  have  been  studied  in 
detail  by  Christodoulidis  (1976).  Probably,  the  variation  in  gravity  due  to  Earth 
tides  and  ocean  loading  should  be  corrected  as  well,  in  order  to  achieve  the 
degree  of  accuracy  required  here.  Another  source  of  systematic  error  is  the 
gravity  net  to  which  the  measurements  are  "tied".  As  explained  later,  system- 
atics  of  more  than  0. 1  mgal  rms  are  undesirable,  so  the  contribution  from  the  net 
should  be  as  small  as  possible.  Master  stations  where  absolute  gravity  is  known 
to,  say,  0.01  mgal  would  be  quite  adequate.  Of  course,  a  constant  bias  due  to  an 
error  in  the  nets'  datum  has  no  effect  on  the  estimated  potential  differences,  and 
can  be  ignored. 

A  further  cause  of  systematic  errors  in  the  gravity  anomalies  are  the 
distortions  in  the  levelling  net  to  which  the  stations  are  tied.  The  influence  of 
such  distortions  on  estimates  of  the  disturbing  potential  have  been  explained 
by  Lelgemann  (1976).  We  are  going  to  study  here  a  way  of  determining  the  anom¬ 
alies  that  minimizes  this  influence. 

Besides  actual  errors  in  levelling,  the  main  reason  for  distortions  in  vertical 
nets  is  the  use  of  tide  gauges  as  benchmarks.  Their  mean  sea  level  marks  are 
supposed  to  be  at  the  same  potential,  and  the  net  is  adjusted  with  this  as  a  constraint. 
In  reality,  the  stationary  sea  surface  topography  already  mentioned  is  present, 
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and  the  potential  differences  among  gauges  are  not  quite  zero.  These  discrepancies 
propagate  as  errors  throughout  the  adjusted  net.  The  larger  the  network,  the 
larger  the  distortions  can  be,  and  also  the  longer  the  distances  over  which  they 
are  correlated.  To  shorten  this  correlation  length  we  shall  break  the  existing 
net  into  smaller  pieces,  and  to  eliminate  the  effect  of  the  sea  surface  we  shall 
not  use  a  net  that  is  adjusted  with  constraints  based  on  tide  gauges.  To  achieve 
this  we  are  going  to  use  the  center  of  each  cap  as  the  levelling  datum  for  all  the 
gravity  stations  inside  that  cap.  The  potential  of  each  station  is  going  to  be 
referred,  accordingly,  to  the  center  point.  Since  the  caps  considered  here  are 
small  ( 5°  and  10°  semi-apertures)  the  levelling  net  for  each  cap  will  consist 
only  of  short  traverses  whose  measurement  errors  can  be  ignored.  If  necessary, 
the  net  inside  each  cap  can  be  adjusted,  to  filter  out  such  errors.  To  use  the 
estimated  potential  of  each  cap  center  in  our  vertical  connections,  we  have  to 
refer  all  of  them  to  some  common  datum.  To  do  this  without  reverting  to  the 
use  of  tide  gauges  for  this  purpose,  we  shall  take  advantage  of  the  accurate 
position  fixes  taken  at  the  cap  centers,  and  find  the  reference  potential  U  at 
each  one  of  them.  If  we  take  the  U  (P.)  for  the  true  potentials,  we  make  a 
mistake  quite  similar  to  that  of  assuming  that  all  tide  gauges  are  on  the  same 
level  surface.  The  errors,  however,  are  proportional  to  the  T(Pi),  the 
quantities  to  be  estimated.  As  shown  below,  this  leads  to  equations  that  can  be 
solved  for  these  "errors",  to  obtain  the  desired  T(Pi)  free  from  biases.  In 
a  way,  it  can  be  said  that  the  centers  of  the  caps  are  the  equivalent  of  tide  gauges 
in  the  adjustment  of  the  World  Vertical  Network. 

A  gravity  anomaly  Ag,  in  terms  of  the  reference  model,  is 


Ag(Q)  =  g(Q)  -  y(Q') 


(3.15) 


where  g(Q)  is  the  acceleration  of  gravity  measured  at  a  point  Q  on  the  Earth's 
surface,  and  y(Q')  is  the  model's  acceleration  at  a  point  Q'  such  that  A.q=Aq< 
and  <pq  =  <pq t,  while  U(Q')  +  «?(Q')  =  V(Q)  +<o(Q)  =  W(Q) .  The  rotation 
potentials  <p(Q)  and  <p(Q')  are  almost  equal,  because  the  difference  ( nj  -  rq» ) 
is  of  the  order  of  3  m  for  a  reference  model  up  to  degree  and  order  20.  For  the 
same  reason,  the  linearized  expression 


Ag(Q) 


-1113)  *y(Q)  T(Q) 

9  r  +  9  r  y(Q) 


(3. 16) 


can  be  regarded  as  almost  exact.  These  approximations  hold  better  for  this  type 
of  model  than  for  the  simpler,  and  traditional,  ellipsoidal  model.  The  gravity 
potential  at  the  gravity  station  Q  is 


W(Q)  =  U(P,) +T(Pt) +AW(P„Q) +^<P0 


where  AW(Pt,Q)  is  the  (levelled)  potential  difference  between  the  cap  center 
Pt  and  Q.  Since  T(Pi)  is  not  known,  we  can  only  measure 


W(Q)  -  T(P,)  =  U(P,)  +  AW(PS,Q)  +  <0(Pt) 


(3. 


where  U(Pt)  is  obtained  from  the  reference  model  and  the  coordinates  of  P, 
(precise  fix).  Therefore,  instead  of  Ag  we  determine 


Ag*  -  g(Q)  -  y(Q") 


(3. 


Q"  is  a  point  with  the  same  p  and  X  as  Q  and  Q' ,  anch  where 


U(Q")  +  <0(Q”)  =  W(Q)  -  T (Pi) 


The  relationships  between  Q,  Q',  and  Q"  are  illustrated  in  Figure  3.  The 
distance  Q"  Q'  is 

CTO5  =  ( U (Qf)  -  U (Q")) y (Q)_1  «  T(P1)y(Q)'1 


and 


So 


y(Q") 


y(Q’)  +  V-(Q')(Q'’Q') 

?r 

y(Q’)  +~(Q)  r(P,)y(Qrx 


Ag*  =  g(Q)  -  y(Q")  -  g(Q)  -  y(Q’)  -  1^2}  T(Pt)  y(Q)"1 

or 


Ag*  =  Ag(Q)  -|*(Q)T(Pt)y(Q)_1  -  Ag(Q)  +”T(p!)  (3. 

The  measured  ag*  can  be  used  in  one  of  two  ways:  (a)  as  a  true  gravity 
anomaly  corrupted  by  a ’TDias"  ^y  T  (Pt)  ; 

5r  y(Q) 


(b)  as  a  function  of  the  gravity  field  on  its  own  right,  defined  by  (3. 19)  as 
dependent  on  both  Ag(Q)  and  T(Pt).  Only  (a)  shall  be  considered  here. 


Figure  3.1.  This  picture  illustrates  the  way  in  which  the  anomalies  Ag 
are  defined  at  the  Earth's  surface,  in  terms  of  the  available 
measurements,  using  the  center  of  the  cap  as  levelling  datum. 
(See  paragraph  (2. 1)  for  notation). 


3.5.  Propagation  of  Position  Errors  through  Gravity  Anomalies 


To  set  up  the  optimal  estimator  (2. 6)  we  have  to  know  the  correlations  among 
the  measurements  and  between  the  estimated  variables  and  the  data.  For  this  we 
use  the  various  correlation  functions  that  are  dependent  on  the  position  of  the  data 
and  the  estimates'  points.  These  functions  are  usually  quite  smooth  horizontally, 
and  an  error  of  up  to  a  few  seconds  of  arc  in  the  geocentric  angle  dioq  is  not  likely 
to  have  any  significant  effect  on  them  and,  thus,  on  the  estimator.  The  same 
functions  are  considerable  more  sensitive  to  a  change  in  the  radial  positions  of 
the  points,  so  errors  in  r  are  more  important  than  those  in  ip.  If  we  have  the 
values  of  Ag  (or  Ag)  at  a  point  Q  of  coordinates  <o,  A,  r,  but  we  assume, 
through  incorrect  position  determination,  that  this  point  is  ^  or  <p  ,  X  ,  r  , 
then 


AAA 

Ag(<p,A  ,  r  ) 


.  SAg  a  * 

Ag(«,A,r)  +  — a  (ns  ,A 
o  r 


A 

r )  er 


?At 


where  eP  =  r  -  r.  The  correction  ^“rB  c r  cannot  be  computed,  as  er  is  unknown, 
so  we  are  confined  to  use  Ag(  Q)  as  if  it  were  Ag(^),  thus  having  a  position  induced 
error  W  fP  in  the  gravity  anomaly.  A  constant  position  error  will  result  in  a 
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constant  bias  in  all  gravity  anomalies,  as  the  variation  of  Tr  over  the  whole  range 
of  r  on  the  Earth’s  surface  (about  ±30  km)  is  very  small.  An  error  in  Ag  will, 
in  turn,  propagate  into  the  estimated  t's 

T  =  £  MAg  +  eAgo  *  T+^eAg,  =  T+0<c  A,  £f,) 

1  =  0  i  1 

If  the  data  arrangement  is  much  the  same  for  each  cap,  then  the  estimator  "weights" 
ft  are  also  much  the  same  in  any  cap,  and  the  resulting  error  in  each  T  from 
a  bias  in  Ag  is  going  to  be  nearly  constant.  These  nearly  equal  errors  in  potential 
will  very  likely  cancel-out  when  potential  differences  are  computed.  For  the  same 
reason,  the  effect  of  an  erroneous  GM  on  the  mean  value  of  the  reference  model's 
gravity  will  not  have  any  appreciable  consequence  on  the  estimated  potential  differ¬ 
ences.  Thus  we  can  have  rather  large  biases  in  GM  and  on  the  set  of  station 
positions  (error  in  the  size  of  the  reference  ellipsoid).  On  the  other  hand,  less 
correlated  errors  in  r  are  not  likely  to  cancel  out.  However,  the  rms  value 
of  at  the  Earth’s  surface  is  approximately  only  8  fj, gals  m”1.  From  the  results 
in  section  4  it  follows  that  errors  up  to  0.1  mgal  in  Ag,  which  are  highly  system- 
matic  inside  a  cap  but  vary  randomly  from  cap  to  cap,  can  be  tolerated.  The  same 
results  show  that  several  mgal  in  errors  totally  uncorrelated  from  station  to  station 
have  very  small  effect  on  the  accuracy  of  T.  From  all  this,  it  can  be  concluded 
that:  a)  errors  of  several  meters  in  r  (at  gravity  stations)  constant  over  the  whole 
Earth,  can  be  tolerated; 

b)  errors  of  several  meters  in  r ,  constant  inside  each  cap,  but  uncorrelated 
from  cap  to  cap,  are  acceptable; 

c)  errors  of  several  meters  in  r,  uncorrelated  from  station  to  station  have 
small  effect. 

The  Global  Positioning  System,  already  mentioned,  is  expected  to  provide  fixes 
accurate  to  10  m  in  each  coordinate  after  six  seconds  of  receiving  radio  signals 
from  four  satellites  visible  at  the  same  time  (Anderle,  1978).  After  15  minutes, 
this  accuracy  should  improve  to  1  m,  plus  a  constant  bias  of  no  more  than  3  m 
due  to  an  error  in  the  adopted  ellipsoid.  Thus,  the  contributions  from  (b)  and 
(c)  above  should  come  to  about  1  m,  guaranteeing  negligible  position  induced 
errors  in  Ag.  By  using  a  GPS  receiver  in  conjunction  with  a  gravimeter,  one 
helicopter  crew  could  collect  all  necessary  information  within  a  5°  cap  (about 
500  stations)  in  less  than  two  months.  Additional  work  would  be  needed  to  estab¬ 
lish  levelling  ties  between  each  station  and  the  prediction  point  at  the  center, 
for  which  already  existing  traverses  could  be  used,  where  available,  jp  fact, 
even  without  GPS  fixes,  the  positions  of  points  like  Q"  in  paragraph  (3.4), 
which  are  obtained  using  levelling  and  the  reference  model,  can  be  used  as  station 
positions.  The  errors  are  going  to  be 

c,q  ~  [  V(Q) -U(Q")]/y  ="  (T  (Pt)  -  T  (Q))/GMa"2 
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If  a  20,20  model  is  used,  the  rms  of  T  will  be  on  the  order  of  3  kgal  m,  so 
*  /2  x  3  kgal  m.  Because  potential  is  a  smooth  function  of  distance,  the 
correlation  c  TT  between  T(Q)  and  T(Pt)  will  reduce  the  rms  of  the  error 

fffr,  =  (Va(T  (Q))*  +  a(T(P))3  -  2  cTT(i/)pQ,  rp.i V,')  v~2 

over  distances  of  the  order  of  one  cap  radius  (500  km),  so  levelling-determined 
positions  might  be  sufficient  for  the  gravity  stations.1 


3.6.  Levelling 


Levelling  is  needed  to  connect  cap  centers  to  benchmarks,  and  cap  centers 
to  gravity  stations  within  their  respective  caps.  Existing  traverses  should  be  used 
wherever  possible,  but  only  unadjusted  values  or  values  adjusted  within  relatively 
small  regions  must  be  used.  Otherwise,  results  could  be  biased  by  the  distortions 
accumulated  in  large  adjusted  nets.  To  keep  errors  small,  the  traverses  should  be 
as  short  as  practicable.  This  goal  is  easily  achieved  for  the  levelling  inside  caps, 
as  distances  from  center  to  rim  range  between  500  and  1000  km  in  the  cases  con¬ 
sidered;  the  longer  traverses  from  centers  to  benchmarks  require  a  careful  planning 
of  the  overall  system.  Caps  should  be  placed,  as  much  as  possible,  within  flat  areas 
with  a  smooth  field  of  gravity  anomalies,  because  estimates  of  T  are  likely  to  be 
better  there  than  where  both  field  and  terrain  change  wildly,  as  suggested  by  the 
results  in  Appendix  A.  Traverses  should  be  levelled  across  regions  where  the 
topography  is  gentle,  to  reduce  their  errors. 

The  simulations  reported  here  have  been  done  assuming  that  all  levelling 
traverses  run  along  arcs  of  maximum  circles,  that  their  errors  are  uncorrelated 
unless  they  overlap,  and  that  the  standard  deviations  of  these  errors  obey  the  simple 
formula 

=  0.  jyXkgalm  (3.20) 

where  £  is  the  length  of  the  traverse  in  thousands  of  kilometers.  This  formula 
represents  a  quality  of  measurement  not  much  better  than  that  of  present  day  first 
order  levelling  (see,  for  instance,  Lelgemann  (1976)). 

An  error  ef|1  will  result  in  an  error  kgal  m  in  U<  P);N  of  about  0.3  erp 

mgal  in  the  value  of  each  A g  ;  and  of  approximately  0. 3  er„  ^fi  kgal  m  in  " 
To.  However,  according  to  Tables  4.3  and  4.4  |  ft  =-0.3  for  S73  and  10°  caps, 
so  <■  To  due  to  erp  is  0(0. 1  crp)  kgal  m  and  can  be  neglected  here  if  fr(<0.5m. 


4.  Computing  Disturbing  Potentials 


Paragraph  (2.2)  explains  the  choice  of  least  squares  collocation  as  the  tech¬ 
nique  for  predicting  T  at  the  cap  centers.  Several  practical  problems  associated 
with  this  method  must  be  considered  before  even  relatively  simple  simulations  can 
be  carried  out.  A  description  of  these  problems  and  their  treatment  is  given  in 
this  section. 


4. 1.  Reducing  the  Dimension  of  the  Data  Covariance  Matrix 


Assuming  that  the  gravity  stations  are  spaced  some  40  km  from  each  other, 
then  their  number  in  each  5°  cap  will  be  close  to  500,  and  to  2000  for  a  10°  cap. 

This  means  that  C«  +  D  is  either  a  500  x  500  or  a  2000  x  2000  matrix,  respectively. 
Being  symmetrical,  it  will  have  up  to  2000000  different  elements  for  a  10°  cap, 
each  one  requiring  one  calculation  of  the  covariance  function  cAi^»  (  P,  Q  )•  This 
can  be  a  very  costly  process  in  terms  of  computer  time,  and  inverting  Ctl  +  D 
can  be  costly  as  well.  This  is  particularly  true  in  the  context  of  a  study  in  which 
calculations  may  have  to  be  repeated  several  times  with  slightly  changed  assump¬ 
tions.  For  these  reasons  an  arrangement  of  the  data  was  chosen  that  gives  the 
C„  +  D  matrix  a  strong  structure.  This,  as  explained  in  Colombo  (1979),  brings 
about  great  savings  in  both  forming  and  inverting  the  matrix.  To  understand  how 
this  is  possible  in  our  case,  consider  the  following  argument.  Imagine  that  the 
gravity  measurements  are  arranged  in  concentric  rows  around  the  center  point, 
which  is  also  the  estimation  point,  and  that  all  stations  are  on  the  same  geocentric 
sphere  as  the  center.  Instead  of  the  individual  anomalies,  suppose  that  we  use 
the  row  sums 

Agt  »  Z  Ag*A 
■  —1 

as  data.  The  number  of  Ag '  s  is  the  same  as  that  of  the  rows,  Nr .  The  covariance 
function  for  the  Ag's  is 

"*! 

C£*£.<M)  s  MfhgiAgj]  =  ^.Z^MfAgt.&gj,}  (4.1) 


(Nt  ,Nj  are  the  numbers  of  stations  in  rows  i,j)  while  that  between  T(  P)  and 
Agi  is 

Ni 

Cr£,(P,i)  =  M(T(P)b£  Ag,.}  =  N;MfTagu] 


as  M(Tagu}  is  constant  for  all  stations  in  the  same  row. 

-21- 


.....  . 


& 


(4.2) 


With  these  two  functions  we  can  set  up  C«  and  C,t ,  while  the  noise  matrix  is 


where  ?ia  is  the  standard  deviation  of  the  mth  measurement  on  the  ith  row.  Cu 
is  now  only  N,  x  Nr  :  for  a  5°  with  40  km  between  rows,  Nr  =  13.  This  implies  a 
reduction  of  more  than  one  order  of  magnitude  in  the  dimension  of  the  data  set  and 
of  several  in  computing  time,  because  setting  up  the  matrix  is  proportional  to  its 
(dimension)  ,  while  inverting  it  requires  a  time  proportional  to  (dimension)3.  More¬ 
over,  as  shown  in  the  reference  given  above,  if  the  points  are  equally  spaced  along 
each  row,  having  the  same  number  in  each  row,  then  estimating  T  from  Ag  or 
from  Ag  gives  the  same  result.  In  other  words:  this  arrangement  largely  improves 
computing  efficiency,  without  changing  the  quality  of  the  result. 

A  distribution  of  data  with  the  same  number  of  points  in  every  row  is  not  a 
good  one,  because  if  the  maximum  separation  (in  the  outer  row)  is  40  km,  then 
the  stations  will  be  very  crowded  at  the  center.  If  points  are  eliminated  to  thin 
out  the  central  region,  the  exact  equivalence  to  collocation  using  Ag  will  be  lost. 

This  will  bring  some  deterioration  in  results,  though  not  a  very  remarkable  one. 

The  actual  accuracy  can  be  found,  as  usual,  using  expression  (2.11)  and  the  matrices 
corresponding  to  £g  's  .  The  unequal  number  of  points  will  result  in  the  outer  rows' 
covariances  being  larger  than  the  central  ones',  and  this  will  probably  worsen  the  con 
dition  number  of  the  matrix.  To  avoid  this  problem,  the  row  sum*  can  be  replaced 
by  row  averages. 

-  1  V 

4g‘ =  -57  .S,18-  <4-4> 


The  covariance  functions  for  these  are 

N, 

CJiS  =  M{Ags  Agj  }  =  jC  M{AglaAg,n}  (4.5) 

Ni 

cr£  =  M[T(P)£g,}  =  E  MfT(P)Ag„}  =  M{T(P)Agtl}  (4.6) 

Nt  »  =  i 

and  ni 

du  =  Vs  E  0?..  D  being  diagonal.  (4.7) 

The  grids  used  in  this  study  were  constructed  according  to  a  simple  pattern  that 
keeps  average  distances  between  stations  at  about  40  km.  There  is  one  station  at 
the  very  center,  six  in  the  first  row,  twelve  in  the  second,  twenty -four  in  the  third; 
their  number  doubling  from  there  on  every  time  the  diameter  of. the  row  doubles  (at 
the  3rd,  6th,  12th,  row,  etc.),  and  staying  constant  otherwise.  Figure  4. 1  shows 
this  scheme  used  on  a  5°  cap.  In  such  a  grid,  the  separation  between  rows  is 


always  almost  40  km,  and  the  constant  separation  of  stations  in  the  same  row  varies, 
from  row  to  row,  from  30  km  to  60  km.  Furthermore,  the  values  of  the  covariance 
between  any  point  in  row  i  and  all  points  in  row  j  2  i  are  repeated  at  all  points  in 
i.  This  greatly  speeds  up  the  creation  of  C«,  as  all  terms  such  as 


"ix 
8  =  1 


MUgi.&g*  1 


in  (4. 5)  are  equal  regardless  of  m.  Finally,  computing  time  can  be  halved  by  taking 
advantage  of  the  fact  that  if  a  radial  line  is  drawn  through  any  point  in  the  figure 
its  covariances  with  all  points  to  the  left  of  this  line  are  the  same  as  with  those  on 
the  right. 


It  is  most  unlikely  that  all  gravity  stations  will  actually  have  the  same  geo¬ 
centric  distances,  and  there  is  no  fundamental  need  for  this,  as  collocation  can  be 
implemented  with  whatever  coordinates  the  stations  may  have,  although  less  effi¬ 
ciently,  as  long  as  they  are  known  with  reasonable  certainty.  However,  if  the  com¬ 
puting  savings  mentioned  above  are  to  be  realized,  the  data  must  first  be  reduced  to  an 
idealgrid  by  collocation  on  a  spherical  surface.  Measurements  in  the  vicinity  of  each  node 
of  the  ideal  grid  can  be  used  to  interpolate  a  value  on  that  node.  Over  a  reasonably 
gentle  terrain,  the  distance  between  it  and  the  sphere  is  not  likely  to  exceed  2  km 
within  a  5°  cap,  and  some  of  the  nodes  are  going  to  be  below,  and  some  above  the 
Earth's  surface.  Assuming  that  five  gravity  stations  were  used,  one  on  the  same 
vertical  as  the  node  but  2  km  above  (below)  it,  and  four  others  forming  cross  with 
the  first  at  the  center  and  20  km  arms,  also  2  km  above  (below)  the  sphere,  the 
accuracy  of  a  value  collocated  on  the  node  is  ±0.8  mgal  if  the  data  has  ±0.5  mgal 
measurements'  white  noise.  This  was  found  using  the  same  covariance  functions 
employed  in  all  the  other  simulations  conducted  during  this  study.  Since  collocation 
is  a  smoothing  process,  the  rms  of  the  estimation  error  is  likely  to  be  due,  by  and 
large,  to  the  high  frequency  components  of  the  data.  As  shown  by  the  simulations 
in  the  next  section,  as  much  as  4  mgal  of  high  frequency  errors  will  have  a  neg¬ 
ligible  effect  on  the  accuracy  of  the  adjusted  vertical  connections. 


4.2,  Estimating  T  from  Ag  instead  of  Ag 


As  already  explained,  the  optimal  estimator 


T (P,)  =  lTd  =  C.z  (C„  +  D)-1  d 


depends  on  the  type  of  data  chosen  d=_z^+£.  Let  Cit  and  Ct:  be  covariance 
matrices  for  gravity  anomalies  Ag  arranged  inside  a  cap  of  center  pB,  but 
suppose  that  the  data  available  consists  in  values  of  Ag  with  the  same  spatial 


-24- 


arrangement.  According  to  (3.19),  if  T  Is  given  in  kgal  m  and  Ag  in  mgal, 


A%  =  Ag  +  k  T(P,) 


(4.8) 


2  *** 

where  k3-  |£- x  106ai  0.3.  Calling  T  the  estimate  of  T(P,)  based  on  Ag,  and 
4*  the  estimate  of  T(F.)  based  on  Ag,  and  using  overbars  to  design  row  averages, 
as  in  the  previous  paragraph,  then 


So 


and 


"r  _____ _________ 

T(P.)  =  T(P,)+7(P.)  =  Eft  [Ag+kT(P.)-k(T(P.))] 

i  =  i 

Nr  Nr 

=  ^fiA^-k^P.)  tu 

T(P,)[1  +kl21f‘]  +  MP.)  =  t5irt  A*g 

Nr  __ 

TCP.)  i  Sift  A* 

5.  \  +  —  —  =  - h  mr 

■  '  ri  .1.  r  f  l  m  il.  f  f 


1=1 


1  =  1  ' 


t<p.) - t,p.i .,(w  - *<*•> ♦  jr;k¥tj 

Consequently, 


N. 


C(P.)  <  €(P.)  if  E  f,  >  o 
1=  1 

Nr 

e(P.)  s  T(P.)  if  E  ft  s  0 


(4.9) 

(4. 10) 

(4. 11-a) 
(4. 11-b) 


1  s  1 


Nr 


If  1^1fi  >0  there  is  a  reduction  in  the  estimate's  error,  according  to  (4. 11-a), 
and  die  opposite  happens  when  I,  <  0  ,  according  to  (4. 11-b). 


4.3.  Accuracies  for  T  ( P.)  Estimated  over  5°  and  10°  Caps 


Tables  4. 1  and  4. 2  show  the  accuracy  of  T  estimated  under  different  con¬ 
ditions,  using  the  theory  explained  in  the  two  preceeding  paragraphs.  To  obtain 
the  theoretical  accuracy,  and  for  the  reasons  given  in  the  previous  paragraph,  the 
square  root  of  the  value  calculated  according  to  (2. 11)  was  corrected  by  the  factor 
( 1  +  k  tSjft)”1.  While  varying  slightly  from  case  to  case,  this  factor  is  always 
close  to  0.9  for  5°  caps,  and  to  0.8  for  10°  caps. 

The  values  of  the  ’’weights"  ft  (i.e.  the  components  of  f)  do  not  change 
greatly,  under  varying  circumstances,  from  the  "typical"  ones  listed  in  Tables 
4.3  and  4.4.  This  is  particularly  true  of  the  largest  "weights",  from  the  center 
point  to  the  10th  ring,  which  remain  nearly  constant. 

In  general  we  can  say  that,  for  a  20, 20  model  and  up  to  several  mgal  rms 
of  white  noise  in  the  gravity  data,  the  accuracy  of  T  estimated  on  a  5°  cap  is 
close  to  0.4  kgal  m,  and  for  a  10°  cap  it  is  near  0. 3  kgal  m. 

The  "imperfect  model"  used  for  the  results  consists  of  the  first  N  degree 
harmonics  of  a  180, 180  model  obtained  by  Rapp  (1978b)  as  a  combination  of  a  world 
data  set  of  l°x  1°  gravity  anomalies  with  GEM-9.  The  standard  deviations  of  the 
coefficients  are  listed  up  to  degree  N  =  30  in  Table  4.5.  The  "2L"  and  "2H" 
covariance  models  are  based  on  formula  (3.8)  and  have  the  following  coefficients 
(jekeli,  1978): 


2L 


2H 


A  =  100  Ofi  =  18.3906  mgal3 
B  =  20  a3  =  658.6132  mgals 
Si  =  .9943667 


sa  =  .9048949 


A  =  140  =  14.0908  mgal3 

B=  10  a3  =  160.6701  mgal2 
sx  =  .9939083 
Sg  =  •  9997595 


where 


.  (Sl 


(f) 


and 


-t 


Covariances  were  computed  with  these  coefficients  and  the  closed  expressions 
also  given  in  the  above  reference.  To  simplify  calculations,  all  measurements 
and  estimations  are  supposed  to  be  made  on  the  same  sphere  of  radius  a  =  6371000 
m.  To  test  the  resulting  accuracies,  some  numerical  experiments  were  conducted, 
as  reported  In  Appendix  A. 
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Table  4. 1. 


Accuracy  of  Estimated  Disturbing  Potential  (kgal  m)  for  5°  Caps  and  "2  L" 
Covariances  (some  values  obtained  using  "2  H"  are  in  brackets). 


Imperfect  Model 

Perfect  Model 

RMS  of  fA% 
in  mgal 

Maximum  Degree, 
Order  in  Model 

0.81 

0.80 

2 

10 

(0.37)  0.41 

- 

4 

20 

(0.36)  0.39 

(0.23)  0.27 

2 

20 

0.38 

0.27 

0 

20 

0.40 

- 

4 

30 

0.37 

0.21 

2 

30 

(Values  estimated  from  Ag's  are  1/0.9  of  those  given  here.) 


Table  4.2. 


Accuracy  of  Estimated  Disturbing  Potential  (kgal  m)  for  10°  Caps 

and  "2  L"  Covariances 


Imperfect  Model 

Perfect  Model 

RMS  of  edfe 
in  mgal 

Maximum  Degree, 
Order  in  Model 

0.29 

- 

4 

20 

0.27 

- 

2 

20 

0.24 

0.19 

0 

20 

i 


Table  4.3. 

Optimal  Row  'Weights"  ft  {kgal  m/mgal)  for  5°  Caps,  Maximum  Degree  and 
Order  in  Model  is  20,  RMS  of  e  A  g  is  2  mgal  (all  values  multiplied  by  0. 1). 


Row  No. 


Imperfect  Model 


"2  L" 


Perfect  Model 


"2  H" 


0  (center) 
1 
2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 


From  these  values  it  is  clear  that  a  constant  error  in  gravity  eg0  (mgal)  will 
result  in  an  estimation  error  (bias) 


«  T0  =  e  go ,  SQfi 


0.3  ego  kgal  m. 


Table  4.4. 

Optimal  Row  "Weights"  fi  (kgal  m/mgal)  for  10°  Caps,  <reAg  =  2  mgal,  imperfect 
Model  to  Degree  and  Order  20,  "2L"  Covariance  Function. 


Row  No.  I  ft  x  0.1  I  Row  No.  !  ftx0.1  !  Row  No.  ftXO.l 


Table  4. 5. 


Accuracies  of  the  Potential  Coefficients  in  the  Imperfect  Model,  by  Degree. 


2899. 

2761. 

2635. 

2522. 

2417. 

2321. 

2232. 

2149. 

2072. 

2001. 


4.4.  Correlation  Among  Estimation  Errors 


The  elements  of  matrix  V  (expression  (2.20))  are  of  the  form 


vklt  =  Mf(€T(PIk)  -  €  T(QJk))3 }  =  (73eT(P1K)+o3fT(QJk)-2MfeT(Plk)cT(QJlt)l(4.12) 


if  the  element  is  on  the  main  diagonal,  or 


Vm  =  M{(C  T(Pt  J  -  €  T(QJB))(c  T(Pi»)  -  (  T(Qj«)) }  (4.13) 

■  MfeT(P,.)«T(P„)}  +  M{fT(Qj.)eT(Qjn)}  -  Mfc  T(P„)  f  T(QJn) }-Mfe  T(Plt)f  T(QJn)  J 

if  it  is  off-diagonal.  To  simplify  matters,  let  us  assume  that  all  caps  are  of  the 
same  size  and  have  the  same  data  distribution  inside.  Then 

*<tVi)  =  * 

where  a  is  the  global  rms  of  the  estimation  errors,  found  in  accordance  to  (2. 11). 
Furthermore, 


The  subscript  k  indicates  the  "caps  pair"  or  observation  equation  number. 


MfcT(ft)  f  T(Qk)}  =  M{(T(R1)-ir^)(T(Qlt)-fTdO)} 

=  c tt ( Rt , Qk)  -  2CrhiltI+ir  C^t  (4.14) 


where  dh  =  [AgltAg  ,...,Agt,...,A®<P]hI  corresponds  to  the  hthcap,  and 

Crhik  =  M{T(Ph)dkr1 

is  a  1  x  Nr  row  vector  (Nr  is  the  number  of  concentric  rows  in  each  cap) 

is  a  Nr  x  Nr  matrix  (<J,  is  the  data  in  the  h  th  cap) . 

The  elements  of  these  matrices  are; 

—  Nl 
Ct  =  M{T(PK)Agl(k)}  =  ^M{T(Ph)  r^Aglr(k)} 

1  r_1 

-  -jj;  r5iM(TtP,,AgL.lki1  (4.15) 


tor  C',», 

(where  Aglr(k)  is  the  value  of  Ag  at  the  rth  gravity  station  on  the  Ith 
ring  centered  at  pk),  and 

"j 

5tr(k)  £  Agj,(h) } 

5  =1 

M{Aglr(k)Agj3(h)]  (4.16) 


'U 


Ni 

*  M{Agt(k) Agj(h)}  =  £  A 

i  J  2» r  =s 
=  NlNTrSx  >4 


for  C^. 

If  the  total  number  of  caps  is  Ne ,  there  are  approximately  (1/4)  Nc2  Nr 2 different 
"ring  covariances"  (the  general  term  in  (4. 16)),  and  each  requires  computing  the 
covariance  function  of  Ag  Nt  x  N4  times.  This  can  tax  the  largest  computing 
budget;  on  the  other  hand  is  only  part  of  the  a  priori  variance-covariance 

matrix  used  for  adjusting  the  potential  difference  between  benchmarks  AW(BMA,BMB). 
Usually,  a  priori  covariances  are  not  needed  to  very  great  accuracy,  because  the 
adjusted  values  are  not  extremely  sensitive  to  them,  or  should  not  be  if  the  procedure 
has  been  designed  properly.  If  approximate  values  are  enough,  then  a  most  efficient 
way  of  obtaining  them,  correct  to  several  significant  figures,  is  to  assume  that  the 
covariance  of  the  "ring  averages"  Ag((h)Ag4(k)  is 
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M{i"*<Mgj<k))  -  «J— L— -  /oSn4g(*,,,)d«7-X;-5Ja"ag(*J,a)d^} 

N 


GSM* 


in^J 


T*-T  £  (2n  +  l)(n-l)afic.Pll<*l)P«rt>dP.(k)  + 

U  n  — a 

N  **  *  s  -1 

+  #£j2n  +  1)(®  -  1)  fin  P»(^l)  P»(0j)  P.(^)  J 


(4.17) 


where  is  the  spherical  distance  between  the  centers  of  the  caps  where  rings 
i  and  j  are  located;  ip  t  and  ip  3  are  the  sizes  (semi-apertures)  of  rings  i  and 
j;  while  N,  6e„,  6„  are  the  same  as  in  expression  (3.4).  Values  computed  as 
above  are  accurate  to  no  less  than  three  significant  figures  fo^  lpt  >  10°,  iplt 
4>]  <  5°  and  N,»*  >  400.  Table  4.6  presents  values  of  Mfe  T(Pk)  eT(Pj) }  for 
different  distances  between  cap  centers:  this  illustrates  the  considerable  indepen¬ 
dence  of  estimates  separated  by  more  than  2000  km.  As  the  data  are  supposed  to 
be  Ag's  rather  than  Ag's,  results  have  been  divided  by  1+0.31!:  fi  ^  0.9 

i  =i 


Table  4. 6. 

Correlation  Between  Disturbing  Potential  Estimates  at  Various  Distances 
(crcAg  =  2  mgal,  5°  Caps,  Imperfect  Model  up  to  Degree  and  Order  20). 


MfeT(P()e 

T(P;)}  (kgal  m)a 

Distance  Pi  Pj 

km 

"2  H" 

"2  L" 

0.130 

0.154 

0 

0.032 

0.033 

1150 

0.027 

0.028 

1300 

0.008 

0.007 

2000 

0.002 

0.002 

2300 

-0.002 

0.002 

2500 

0.001 

0.002 

14000 

-0.002 

-0.002 

17500 

5.  The  Accuracy  of  the  Adjusted  Vertical  Connection 

This  section  presents  the  main  results  of  this  study;  the  theoretical  accuracy 
of  a  World  Vertical  Network  constructed  along  the  lines  already  discussed.  It  also 
contains  the  theory  of  an  optimal  estimator  for  potential  differences  between  centers 
of  pairs  of  identical  caps,  and  other  ideas  regarding  possible  uses  of  satellite  and 
terrestrial  data  for  setting  up  and  strengthening  levelling  nets. 


5.1.  Transoceanic  Connections  Using  Several  5°  Caps 


Two  cases  have  been  considered:  in  the  first,  five  5°  caps  were  placed 
in  North  America  (4  in  the  U.S.  and  1  in  Canada)  and  another  four  in  Australia. 
The  cap  centers  have  been  chosen  to  avoid  overlaps  and  to  ensure  that  almost 
all  the  area  covered  by  the  caps  is  land.  The  "benchmarks",  two  points  whose 
coordinates  are  given  below,  are:  BMA  in  Electra,  Texas,  and  BMB  in  Wiluna, 
Western  Australia.  Cap  centers  in  the  same  land  mass  are  supposed  to  be  joined 
to  the  corresponding  benchmark  by  levelling  traverses,  as  in  Figure  2. 1.  To 
assign  a  standard  deviation  to  the  error  of  each  traverse,  formula  (3.20)  has 
been  used,  the  length  of  the  traverse  being  that  of  the  maximum  circle  joining 
its  endpoints.  All  data  is  supposed  to  have  been  measured  (or  reduced)  to  the 
same  sphere  of  radius  a  =  6371000  m  to  simplify  calculations.  Tables  5. 1  and 
5.2  show  the  accuracies  obtained  using  formula  (2. 19;  with  the  "a  priori"  matrix 
V  corresponding  to  various  combinations  of  covariance  function  and  reference 
model.  The  following  is  the  list  of  cap  centers,  including  their  latitudes  (the 
benchmarks  are  also  cap  centers).  The  difference  between  the  North  America/ 
Australia  and  the  U.S. /Australia  connections  is  the  fact  that  cap  number  9  in 
the  former  has  been  replaced  by  cap  9*  in  the  latter. 


Cap  No. 

Latitude 

Longitude 

Location 

State 

1  (BMB) 

-27.5° 

120.0° 

WILUNA 

W.  Australia 

2 

-20.0° 

130.0° 

Tana  mi 

N.  Territory 

3 

-22.5° 

142.0° 

Middleton 

Queensland 

4 

-32.5° 

145.0° 

Cobar 

N.  S.  Wales 

5 

36.0° 

-85.5° 

Sparta 

Tennessee 

6  (BMA) 

34.5° 

-99.0° 

ELECTRA 

Texas 

7 

44.0° 

-99.0° 

Wessington 

South  Dakota 

8 

37.0° 

-112.5° 

Fredonia 

A  rizona 

9 

56.5° 

-112.0° 

Me.  Murray 

Alberta 

9* 

65.0° 

-150.0° 

Manley 

Alaska 

Overall,  the  accuracies  listed  in  Tables  5.1  and  5.2  can  be  separated  into  two 
groups:  those  based  on  the  use  of  an  imperfect  reference  model,  and  those 
obtained  assuming  a  perfect  model.  Results  within  each  group  are  much  the 
same:  approximately  0.  3  kgal  m  accuracy  with  an  imperfect  model,  about  0. 2 
kgal  m  with  a  perfect  one.  Changes  of  100  in  the  accuracies  of  levelling  or 
gravity  measurements  caused  less  than  10%  variation  in  a  AW  (the  effect  of  lev¬ 
elling  accuracy  is  shown  in  Table  5.2);  while  a  change  in  covariance  model  had 
only  a  slight  effect  on  AW ,  as  suggested  by  the  components  of  the  respective 
pseudolnverse  vectors  v  listed  in  Table  5.3  (where  vr£  =  AW,  see  paragraph 
(2.4)),  and  as  shown  in  Table  5.1. 
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In  summary:  the  accuracy  of  the  reference  model  is  the  most  important  of  the 
various  sources  of  error  included  in  this  study.  Improvements  in  this  model 
are  likely  to  have  a  large  effect  on  the  quality  of  the  resulting  vertical  network. 

As  the  error  correlations  between  caps,  shown  in  Table  4.  6,  are  very 
small  compared  to  the  autocorrelations  (cap  errors),  we  can  ignore  them  in 
a  first  approximation,  regarding  all  errors  contributing  to  matrix  V  as 
uncorrelated.  The  standard  deviation  of  the  adjusted  iW  can  then  be  guessed 
using  the  following  formula,  instead  of  the  "exact"  expression  (2. 19): 


<7  AW  (BMA,  BMB)  ( 2  ( C^CT  +  0.01  4,»x)  +  y2Oae  A  r)^/  N e‘ 

where  Ne  is  the  number  of  independent  equations  (paragraph  (2. 4)),  £»,x  is 
the  length  of  the  longest  traverse  (  <  5000  km),  and  ore  Ar  is  the  error  in  rela¬ 
tive  geocentric  position  between  points.  Assuming:  eight  equations,  as  in  the 
two  examples;  a*t  =  2  mgal;  and  the  imperfect  model,  then 


and 


crAW(BMA,BMB)  ^  (0  .  40  +  0.96  o2eAr)V/8  for  "2L" 
<TAW(BMA,  BMB)  =-  (0.36+0. 96aseAr)^//8  for  "2  H" 


Assuming,  for  instance,  that  (7€Ar  %  0.5  m,  the  corresponding  accuracies  with 
these  simplified  expressions  would  be  0.28  kgal  m  and  0.27  kgal  m,  respectively. 
Now  ±0.5  m  in  relative  geocentric  position  is  within  the  reach  of  present-day 
Doppler  satellite  techniques,  hi  any  case,  all  the  results  given  here  are  well 
below  the  ±1.5  kgal  m  theoretical  uncertainty  for  tidal  gauges  and  spirit  levelling 
alone. 


Table  5. 1. 


Accuracy  <7AW  of  the  North  America-Australia  Vertical  Connection  (in  kgal  m) 
o,eAr  =  0.15  m,  oeAg=2mgal,  ere AW£  =  0.  l/l  kgal  m 
________________________ 


Imperfect  Model  (N  =  20) 


Perfect  Model  (N  =  20) 


Table  5.2. 


Accuracy  aAW  of  the  U.S.  A. -Australia  Vertical  Connection  (in  kgal  m) 
2  L  covariance,  aeAg=2mgal,  acAr  =  0.15  m 


Imperfect  Model  (N  =  20)  Perfect  Model  (N  =  20)  crcAW 


Table  5.3. 

Components  of  v=  (aT  V-1a)_1aT  V"1  (Dimensionless) 

North  America-Australia  Connection 
crfAg  =  2  mgal,  creAr  =  0.15  m,  crcAW^  =  0. 1/ikgal  m 


"Observation  Equations” 

Potential  Difference 
between  caps  number 

Eqn.  No. 

1-5 

1 

1-6 

2 

2-6 

3 

2-7 

4 

3-7 

5 

3-8 

6 

4-8 

7 

4-9 

8 

Imperfect  Model  (N  =  20)  Perfect  Model  (N  =  20 


•rot"  MO  U" 


!»9  T  If  tTO  TTM 


.207 

.204 

.048 

.050 

.204 

:  .212 

.098 

.108 

.104 

.085 

.137 

.155 

.072 

.044 

.129 

.140 

timal  Estimator  for  the  Potential  Difference  Between  Two  Caps 


In  addition  to  the  more  general  configuration  studied  in  the  preceeding  para¬ 
graphs,  we  shall  consider  a  ’’minimal  estimator”  where  only  two  caps  are  involved, 
each  centered  at  a  benchmark.  As  in  paragraph  5. 1,  we  shall  restrict  the  esti¬ 
mator  by  assuming  that  the  data  Ag  has  been  converted  to  ring  averages  5g  and, 
furthermore,  that  both  caps  are  identical  in  size  and  data  arrangement.  This  limits 
somewhat  the  power  of  the  estimator,  in  particular  the  use  of  ring  averages  makes 
it  suboptimal  by  comparison  to  a  "full”  estimator  based  on  point  data.  On  the  other 
hand,  these  constraints  greatly  expedite  creation  of  the  estimator,  as  in 
the  case  considered  before. 


'..I 


The  objective  of  the  estimator 


A  ft  NP 

AT(Pi,Ps)  =  -  fid,  =  l  fu  Agi(l)  -  Z  fjaA"gj(2)  (5.1) 


is  to  minimize  the  global  mean  square  of  the  prediction  error: 


M{fA"T(P1,  P3  )a  }  =  M{(AT(PlfPa) -fjda))8} 


Because  both  caps  are  identical,  and  the  average  is  isotropic  (function  of  spherical 
distance  and  radial  distances  only)  then,  if  both  caps  are  on  the  same  sphere,  the 
optimal  weights  for  each  are  the  same; 

f)  =  _f  2  =  _f  or  f tl  =  fl  2  =  fl  i  =  1.  2,  ...»  Nr 

Accordingly, 


*MMT(Pi,Pa)^  =  ^M{(T(P1)-T(^)  -  fi(Agi(l)  -  Agi(2) )  )a] 

«  iM{(T(Pt)-T(Ii)  -  f\&-djJ))<T(Pl)-T(Pa)-(4-dfl)T.f)  1 
=  crr(0)  -Ctt(^i%)  +  2J[ T(Cf]jt2-Crl4l )  z2+D)_f 


where  cz^  is  the  same  as  Ca^  in  (4. 14),  and  D,  Ct^  =  Ct^  areas 

in  (2.7).  Thus, 


i  ~  M  f  CAT)  =  Ctu  -Cr^.  +(CSl,1-CH.,+  D)f  =  0 

Cl 


_f  r~ C 1 2+  D)  (Cr^ "CTjdj) 


since  the  data  consists  in  values  of  Ag  rather  than  Ag,  a  "correction  factor" 
should  be  used,  as  before 


it  -  L'(4gi  •  i&>  *  x+o.aj.t,  <5-5 

where  , . . .  ,Z"gNr  ]T.  Using  the  type  of  data  grid  shown  in  Figure  4. 1, 


extended  to  a  10°  cap,  so  that  the  average  separation  of  the  stations  is  40  km, 
assuming  the  same  reference  gravity  fields  and  covariance  models  (2  L  and  2  H) 
used  before,  0,1m  relative1  position  accuracy  (cap  centers)  and  2  mgal  accuracy 
for  the  gravity  anomalies,  the  following  global  rms  of  the  estimation  errors  were 
calculated  for  various  separations  $4  between  the  caps: 

Table  5.4. 

RMS  of  Error  of  Optimally  Estimated  Potential  Difference  Between  the  Centers  of 
Two  10°  Caps  1^4°  Apart.  (Imperfect  model  to  degree  and  order  20,  <7eA%  =  2  mgal, 
"2  L"  covariance  function. ) 


(J£iW 

(kgal  m) 

<°> 

0.35 

10° 

0.40 

20c 

0.42 

30° 

(RMS  fluctuates  between  0.41  and  0.42  kgal  m  for  30°  <  180°.) 

If  the  position  errors  at  each  cap  center  are  uncorrelated,  and 
rms,  then  a  (  AW  should  be  corrected  as  follows;  cr  e  AW'  = 

It  is  interesting  to  compare  Table  5.4  to  Table  4.2:  for  itid  > 
listed  above  clearly  approaches  Jz  MfeT3]  where  /MffT2)  is  the  "single  cap" 
accuracy  listed  in  Table  4.  2.  This  also  agrees  with  the  increasing  independence  of  es¬ 
timates  of  T  separated  by  more  than  a  few  thousand  kilometers,  pointed  out  in 
paragraph  4. 4.  Similarly,  the  values  of  the  optimal  "weights"  ft  converge  to 
those  for  a  single  10°  cap,  with  increasing  ^4. 


5.3.  Height  Differences  Between  Inaccessible  Points 


The  technique  described  in  this  report,  in  essence,  uses  accurate  position 
fixes  and  a  good  gravity  field  model  to  obtain  an  estimate  of  the  potential  difference 
between  points,  estimate  that  is  then  refined  using  data  from  the  neighborhood  of 
each  point  (gravity  anomalies  and  levelling)  to  add  high  frequency  information  not 
contained  in  the  model.  So  far  we  have  assumed  a  good  local  coverage,  easily 
accessible  gravity  stations  and,  generally,  a  most  cooperative  disposition  from 
both  Nature  and  men  towards  our  project.  If  either,  or  both,  were  lacking,  we 
would  be  left  with  the  field  model,  plus  some  data  in  the  perifery  of  the  region  of 
interest  to  refine  the  former,  and  perhaps  not  even  a  high  accuracy  relative 
position  fix  at  each  cap  center.  Through  collocation,  external  data  can  be  incor- 
porated  into  the  adjustment,  though  not  as  efficiently  as  in  the  scheme  discussed 
1  The  relevant  error  here  is  that  in  the  measured  difference  of  geocentric  distances. 
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in  paragraph  4. 1.  As  the  data  is  far  from  the  point  of  interest,  it  is  important  to 
use  information  rich  in  long  wavelengths  signal.  One  possibility,  if  the  inacces¬ 
sible  area  is  not  extremely  far  from  the  sea,  is  the  mean  sea  surface  derived  from 
satellite  altimetry,  regarded  as  a  quasi-geoid  with  1,5  m  global  rms  of  ”noise" 
due  to  dynamic  effects.  Altimetry  could  be  used  on  land  as  well,  to  provide  fixes 
in  radial  position  for  the  "cap  centers".  The  surface  of  an  inland  sea  or  lake 
would  be  ideal  as  a  target;  other  areas  where  consistent  ellipsoidal  heights  could 
be  obtained  (i.  e.,  independent  of  seasonal  effects),  such  as  salt-flats  in  desert 
areas,  etc.,  could  be  used  if  systematic  errors  due  to  surface  reflectivity  were 
sufficiently  understood.  With  very  accurate  gravity  field  models  and  continuous 
tracking  of  the  altimeter  satellite  by  another  craft  in  a  higher  orbit  it  should  be 
possible,  eventually,  to  obtain  computed  orbits  good  to  0.1  m  (rms),  which,  added 
to  another  0.1m  (rms)  error  for  the  altimeter  itself,  would  amount  to  near  0. 15 
m  (rms)  error  in  the  ellipsoidal  height  fixes,  or  about  0.2  m  relative  geocentric 
height  -^curacy  between  benchmarks  which  would  propagate  as  a  0.2  kgal  m  error 
in  potential  difference.  Gravity  field  models  have  been  improving  at  a  fast  rate 
in  recent  years;  new  and  better  tracking  systems  are  being  developed;  fall 
coverage  of  altimetry  data  over  the  oceans  is  now  available  from  the  Geos-3  and 
Seasat-1  spacecrafts.  These  are  three  encouraging  signs  that,  in  the  coming 
decade,  there  will  be  enough  information  to  model  the  field  up  to  degree  180  with 
a  global  residual  rms  of  about  1  kgal  m  for  the  disturbing  potential.  Then,  vertical 
connections  good  to  at  least  /2  x  1  +  (0.2)*  =“1.5  kgal  m  should  become  feasible. 
This  is  the  same  as  the  theoretical  global  accuracy  of  a  system  based  on  tide 
gauges;  it  should  be  much  the  same  as  having  "tide  gauges"  inland. 


5.4.  Some  Questions  Regarding  Accuracy  Estimates 


The  accuracies  listed  in  the  various  tables  of  this  section  depend,  mostly, 
on  the  applicability  of  the  theory  of  collocation  to  the  real  world.  In  particular, 
there  are  some  aspects  of  collocation  open  to  criticism  that  deserve  a  mention 
here.  The  first  one  is  the  unavoidable  use  of  approximate  covariance  functions,  as 
the  true  ones  cannot  be  known  exactly  from  finite  amounts  of  data,  basically  because 
the  expansion  of  the  "true"  gravity  field  in  harmonic  polynomials  is  infinite.  Laurit- 
zen  (1973)  has  established  the  impossibility  of  obtaining  the  covariance  of  a  random 
process  on  a  sphere  even  from  a  complete  data  coverage,  but  it  is  not  very  clear 
why  the  gravity  field  should  be  treated  as  a  random  process  in  the  first  place.  Be 
as  it  may,  the  available  information  is  always  going  to  be  incomplete  and  inaccurate. 
Otherwise,  there  wouldn't  be  much  point  in  using  collocation,  or  any  other  form  of 
interpolation  and  filtering,  as  there  would  be  little  to  learn  from  it.  This  is  why 
two  different  empirical  covariance  functions  were  used,  "2  L"  and  "2  H",  to  find 
out  how  sensitive  the  adjusted  potential  differences  were  to  the  choice  of  function. 

The  results,  as  shown  in  the  various  tables,  were  much  the  same  with  both. 


Second,  there  is  the  question  of  how  suitable  is  the  spherical  harmonic 
representation  of  the  gravity  field  on  which  a  good  deal  of  the  theory  put  forward 
here  (as  much  of  modern  geodesy)  rests.  From  the  work  of  Petroskaya  (1977) 
and  Sjoberg  (1978),  who  have  proposed  mathematical  formulations  for  the  field 
both  inside  and  outside  the  Earth's  surface,  we  learn  that  such  expansions  are 
not  necessarily  sums  of  solid  spherical  harmonics.  On  the  other  hand,  it  is 
well  known  that  the  gravity  field  of  a  homogeneous  ellipsoid  can  be  expressed  in 
terms  of  such  harmonics  down  to  a  sphere  completely  buried  in  the  ellipsoid. 

The  Earth  being  primarily  ellipsoidal,  we  may  expect  it  to  have  a  field  that  does 
not  behave  too  differently  from  that  of  the  ellipsoid.  However,  there  is  no  reason 
to  expect  that  the  harmonic  series  that  describes  the  field  exactly  outside  the 
bounding  sphere  does  so  also  inside  it  and  down  to  the  Earth's  surface.  That  it 
does  not  diverge  too  strongly  is  shown  by  the  fact  that  low  frequency  models  de¬ 
rived  from  satellites  provide  a  fit  to  surface  data  that  improves  as  more  terms 
are  used.  That  the  model  is  not  too  inadequate  follows  from  the  usefulness  of 
formulas  such  as  Vening  Meinesz'  and  Stokes',  which  are  based  on  the  assumption 
that  the  field  can  be  expanded  in  solid  spherical  harmonics.  But  all  this  supporting 
evidence  shows  is  that  these  ideas  "work"  enough  to  provide  about  one  meter  ac¬ 
curacy  in  computed  undulations,  seconds  of  arc  in  deflections  of  vertical,  a  few 
milligals  in  interpolated  gravity.  No  evidence  is  available  to  suggest  that  they 
also  "work"  at  the  level  of  accuracy  (approximately  0.3  to  0.5  m)  expected  of  them 
here.  Going  back  to  the  homogeneous  ellipsoid,  it  is  sufficient  to  add  a  most  minute 
inhomogeneity,  in  the  form  of  a  tiny  material  sphere  at  any  distance  R  from  the 
origin,  for  the  series  of  the  composite  body  to  fail  to  converge  inside  the  sphere 
of  radius  R.  The  mass  of  the  sphere  (and  therefore  the  difference  between  the 
pure  ellipsoidal  and  the  composite  field)  can  be  made  as  small  as  desired  without 
the  series  ever  converging  again.  This  makes  clear  that,  at  least  in  some  cases, 
a  gravity  field  with  a  harmonic  series  that  does  not  converge  down  to  an  internal 
sphere  can  be  only  slightly  different  (except  at  some  isolated  points)  from  one  with 
a  series  that  does.  Krarup  (1969)  brought  attention  to  this  fact,  and  enquired 
whether  this  might  not  be  always  the  case.  He  concluded  that  it  is,  furnishing 
proof  of  what  he  called  a  "Runge-type  theorem",  because  of  similar  theorems  for 
elliptical  differential  equations.  Krarup's  thesis  is  that  any  harmonic  field  can  be 
approximated  uniformly,  together  with  all  its  derivatives,  by  sequences  of  series  of 
spherical  harmonics  that  converge  to  an  arbitrarily  small  sphere  centered  at  the 
origin,  completely  inside  a  surface  (terrain)  that  separates  the  region  where  the 
field  is  harmonic  from  that  where  it  is  not.  There  are  a  few  restrictions  on  the 
nature  of  this  surface,  but  they  are  loose  enough  to  ensure  an  adequate  fit  to  the 
real  topography.  The  uniform  convergence  takes  place  only  down  to  that  surface. 
Inside,  the  various  approximations  may  differ  increasingly  from  each  other,  and 
have  no  limit  function.  So,  a  spherical  harmonic  approximation  to  the  exterior 
field  can  be  quite  irregular  and  "wild"  in  the  interior,  and  this  might  be  a  cause 
for  some  concern  here.  The  global  rms,  or  accuracy  obtained  from  (2.11),  cor¬ 
responds  to  the  average  of  the  square  of  the  errors  of  all  possible  predictions  made 
on  the  same  sphere  where  the  actual  estimation  point  is  situated.  Such  sphere  is 
always  partly  inside  the  solid  Earth,  because  of  the  equatorial  bulge.  If  we  imagine 
the  covariance  functions  that  we  are  using  as  corresponding  to  some  spherical  har¬ 
monic  expansion  that  fits  closely  the  field  outside  the  terrain,  they  may  also  cor- 
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respond  to  a  function  that  has  a  markedly  different  character  on  that  part  of  the 
sphere  that  is  buried  from  that  which  is  out  in  the  open.  In  other  words,  the 
spherical  harmonic  approximations  may  not  be  sufficiently  "stationary"  for  a 
meaningful  application  of  collocation.  It  may  help  to  understand  this  problem 
having  some  method  to  generate  the  sequences  of  convergent  series  Krarup  has 
spoken  of,  in  order  to  visualize  their  behavior.  Such  method,  to  this  author's 
knowledge,  has  not  been  proposed  yet. 

The  real  question  here  is  whether  collocation,  and  the  use  of  spherical 
harmonics  theory,  may  not  be  a  source  of  bias  in  the  results.  In  practice  there 
may  be  many  sources  of  bias  not  treated  here,  such  as  systematic  errors  in 
gravity  measurements,  position  fixes,  etc.  The  only  definitive  way  of  knowing 
if  such  factors  can  be  truly  significant  is  to  subject  the  whole  idea  to  experimen¬ 
tation.  If  enough  determinations  of  potential  differences  are  carried  out  using 
this  technique,  at  many  places  around  the  world,  between  points  already  connected 
by  levelling  traverses,  then  any  large  systematic  discrepancies  should  become 
apparent.  If  they  do  not  show  up,  then  the  experiments  will  provide  supporting 
evidence  for  the  use  of  the  idea  elsewhere. 


6.  Conclusions 

This  report  has  dealt  with  the  concept  of  World  Vertical  Network:  a  set  of 
benchmarks  thousands  of  kilometers  apart,  the  geocentric  coordinates  of,  and  the 
potential  differences  between  which,  are  accurately  known,  so  they  can  be  used  to 
connect  separate  levelling  nets  in  a  world-wide  system. 

The  results  of  section  5  suggest  that,  according  to  theory,  it  might  be  possible 
to  set  up  this  network  to  a  significantly  better  accuracy  than  that  provided  by  tide 
gauges  and  spirit  levelling  alone.  This  could  be  done  by  using  satellite  and  terres¬ 
trial  data  together,  employing  least  squares  collocation  as  the  main  mathematical 
tool  for  combining  them. 

Two  cautionary  notes  are  in  order:  first,  there  may  be  causes  of  error  not 
considered  in  this  study  that  turn  out  to  be  of  practical  significance*  only 
actual  experience  can  have  the  final  word  on  this.  Second,  the  estimates  of  T, 
crucial  to  the  ideas  in  this  work,  are  assumed  made  from  values  of  Ag  measured, 
or  interpolated,  on  spherical  surfaces.  Collocation  can  be  used,  in  principle,  on 
more  general  (terrain-like)  surfaces,  but  the  accuracy  is  not  exactly  the  same 
as  with  the  spherical  arrangement  of  data  points.  "Common  sense"  suggests  that, 
if  the  terrain  is  gentle,  departing  smoothly  from  the  spherical  shape,  the  accuracy 
should  be  much  the  same,  but  no  evidence  for  this  is  given  here.  Numerical 
studies  are  far  more  difficult  when  the  data  is  not  on  a  simple  spherical  arrange¬ 
ment,  and  probably  too  expensive  for  an  ordinary  research  project.  Furthermore, 
whether  on  the  terrain  or  on  a  partially  buried  sphere,  the  statistical  relevance  of 
the  accuracies  derived  from  collocation  can  be  questioned  on  the  grounds  given 
at  the  end  of  paragraph  (5.4).  If  this  is  a  real  problem,  it  is  in- 
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he  rent  to  all  applications  of  collocation  theory  to  estimation  of  gravitational 
field  variables  inside  the  Earth's  bounding  sphere,  and  not  just  to  the  present 
ideas.  This  matter  seems  to  bu_ve  received  little  or  no  attention  from  workers 
in  this  area. 

Depending  primarily  on  the  quality  of  the  reference  gravity  field  model,  the 
accuracy  of  the  vertical  connections  in  the  network  could  be  between  0. 2  and  0.3 
kgal  m,  using  the  data  arrangement  described  in  sections  2  and  5.  Such  configu¬ 
ration  has  been  chosen,  mostly,  to  simplify  this  study,  and  is  by  no  means  the 
only  possible  one.  With  more  data,  including  satellite  altimetry,  we  could  expect 
even  better  results,  so  those  given  here  are  to  be  regarded  as  upper  limits  to  the 
ultimate  quality  of  a  global  network. 

The  North  A  meric  a- Australia  connection  studied  in  section  5  requires 
relatively  few  data  (about  4000  point  gravity  anomalies  to  2  mgal  accuracy,  8 
accurate  position  fixes  and  a  number  of  levelling  traverses,  plus  a  reference 
model  to  degree  and  order  20)  and  the  quality  of  the  measurements  are  almost 
entirely  within  present  day  limits.  The  exception,  the  position  fixes,  can  be  ex¬ 
pected  to  become  feasible  within  the  coming  decade.  Much  of  this  data  can  be 
obtained  and  used  for  other  purposes,  such  as  the  study  of  polar  motion  and  Earth 
rotation  in  the  case  of  the  accurate  point  coordinates.  In  this  way,  by  sharing  with 
other  scientific  enterprises,  the  establishment  of  the  World  Vertical  Network  could 
be  made  both  cheaper  and  an  integral  part  of  the  creation  of  a  World  Geodetic 
System  for  positions  and  heights. 
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Appendix  A 

Verifying  the  Correctness  of  the  Accuracies  Given  in  Section  4 


The  covariance  matrix  for  ring-averaged  gravity  anomalies  is  very  poorly 
conditioned,  with  the  arrangement  of  Figure  4.1,  because  of  the  closeness  of  the 
measurements.  Evidence  of  this  was  observed  while  trying  to  use  tables  with 
equispaced  entries  and  linear  interpolation,  as  an  inexpensive  way  of  computing 
the  covariance  function  values  needed  to  set  up  the  matrix.  With  spacings  as 
small  as  1  km,  the  small  perturbation  due  to  interpolation  errors  was  enough 
to  produce  a  matrix  with  negative  eigenvalues,  which  the  true  covariance  matrix 
can  never  have.  On  the  light  of  this  experience,  it  appeared  reasonable  to  ques¬ 
tion  whether  any  results  dependent  on  such  a  ticklish  matrix  could  be  regarded  as 
meaningful,  including  the  accuracy  estimates  of  Tables  5.1  and  5.2.  To  clarify 
this  matter,  a  very  high  degree  field  model  (up  to  n  =  1000)  consisting  only  of 
zonal  terms  was  used 


~  rvr  1000 

T (&)  =  £  c„  P„(0)  (A.l) 

n  3=  a 

with  the  cn  selected  so  the  disturbing  potential  T  would  have  the  same  spectrum/ 
covariance  function  as  the  one  assumed  for  the  Earth.  Using  the  appropriate  ex¬ 
pansion  for  the  gravity  anomaly  of  this  model,  Ag  was  calculated  at  every  point 
the  grid  of  Figure  4.1,  and  then  the  ogtimal  "weights"  from  Table  4.4  were  used 
obtain  the  estimated  The  error  er  =  ^  -  T  was  then  calculated  using  once 
more  the  model,  and  the  whole  operation  was  repeated  at  5°  intervals,  from  pole 
to  pole.  The  mean  value  of  the  squared  error  is,  approximately 

9,=n 

M[Tr8}  =“  2*ra  £  sin  Sj  MS,  )a/(2rra  £  sinSt)  (A. 2) 

01=0  01=0 

-  ^  A  ^ 

where  advantage  is  taken  of  the  fact  that  the  error  fT  as  well  as  T,  T?  and  Ag 

have  expansions  consisting  of  zonal  terms  only.  Assuming  a  perfect  reference 
model  up  to  degree  and  order  20,  2  mgal  (rms)  error  in  the  gravity  anomalies, 
the  global  rms  of  the  prediction  error  (accuracy),  estimated  using  (A. 2)  was  0. 3 
kgal  m.  This  is  about  30%  below  the  value  given  in  Table  4. 1,  but  not  unreasonably 
so*  the  sampling  near  the  North  pole,  at  5°  intervals,  is  probably  too  coarse  to 
accurately  cover  the  fast  changes  in  the  field  there.  The  largest  errors  also 
occur  near  that  pole,  so  a  shorter  interval  is  likely  to  pick  up  more  of  them,  in¬ 
creasing  the  right  hand  side  of  (A.  2). 

Table  A.  1  shows  T  and  the  errors  Cr  ,  as  functions  of  latitude 
(10°  intervals).  The  largest  errors,  as  already  mentioned,  coincide  with  the  wild 
"spike"  in  the  disturbing  potential  near  the  North  pole.  The  error  outside  a  30° 
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polar  cap  is  almost  everywhere  less  than  the  global  rms.  This  suggests  that  the 
worst  predictions  can  occur  where  the  field  has  large  and  fast  variations.  The 
percentage  of  error  in  the  predicted  value,  however,  is  about  10%  at  the  pole,  so 
the  relative  goodness  of  the  prediction  is  not  particularly  bad  there.  But  we  are 
interested  in  the  absolute  value  of  the  error,  so  the  pole  is  clearly  a  bad  place, 
in  this  "artificial  planet",  for  making  estimates.  Table  A.  1  shows  that  the  be¬ 
havior  of  Ag  closely  resembles  that  of  ^  .  In  this  particular  case,  the  predic¬ 
tion  of  is  bad  where  that  of  2g  is  bad,1  and  vice  versa,  and  this  could  be 
used  as  a  practical  criterion  for  selecting  the  locations  of  estimation  points. 
Predicting  Ag  from  gravimetry  in  a  "candidate"  region,  at  points  in  that  region 
where  Ag  is  already  known  from  accurate  measurements,  we  could  compare 
the  actual  error  in  these  estimates  to  the  theoretical  accuracy  of  expression 
(2.11).  If  the  actual  errors  are  close  to  their  theoretical  rms,  the  region  is 
acceptable,  and  we  can  proceed  to  set  up  a  cap  like  that  one  in  Figure  4. 1;  if  the 
errors  are  consistently  larger,  the  region  should  be  rejected. 

Other  points  to  consider  when  selecting  a  place  for  a  cap  are:  lack  of  sig¬ 
nificant  bodies  of  free  water,  because  gravity  cannot  be  measured  as  accurately 
on  water  as  on  land;  smooth  topography  (with  most  places  easily  accessible)  to 
permit  accurate  levelling  and  good  coverage  with  gravity  stations. 


Ag  is  estimated  here  using  a  "4  points'  cross"  pattern  centered  at  the  prediction 
point.  The  "arms"  are  20  km  wide. 
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Table  A.  1 


Prediction  Errors  in  T  and  Ag  in  a  Zonal  Field.  Imperfect  Reference 

Field  to  Degree  and  Order  20,  5°  Cap,  ae  Ag  (  2  mgal  (  T) 

l  0.5  mgal  ( Ag) 


<0 

<°) 

A 

(T 

kgal  m 

mgal 

T 

kgal  m 

Ag 

mgal 

90 

80 

-72.88 

.53 

3195. 39 
-1.24 

697.91 

-3.60 

26095.57 

-14.36 

70 

.06 

.20 

-  .50 

-  2.14 

60 

-.17 

.20 

.77 

4.15 

50 

-.31 

.37 

-1.26 

-  3.20 

40 

.04 

-  .28 

.78 

2.41 

30 

-.  13 

-  .07 

-  .60 

-  1.81 

20 

-.01 

-  .04 

.13 

.45 

10 

-.02 

.08 

.01 

.60 

0 

-.06 

.14 

-  .33 

-  .89 

-10 

.06 

-  .12 

.41 

1.17 

-20 

-.02 

-  .03 

-  .34 

-  1.34 

-30 

.07 

-  .07 

.31 

.80 

-40 

.04 

.96 

-  .03 

-  .15 

-50 

.05 

.10 

-  .03 

-  .15 

-60 

.07 

-  .09 

.36 

.97 

-70 

-.12 

-  .04 

-  .58 

-  1.85 

-80 

-.11 

-  .20 

.46 

1.92 

-90 

-.69 

3.63 

-2.25 

-  .65 

This  Appendix  contains  listings  of  parts  of  the  software  developed  for  this 
project,  and  sample  output.  First  there  is  a  listing  of  "RINGS",  a  program  that 
obtains  the  optimal  components,  or  weights,  of  the  estimator  vector  _f  (expression 
(2.6)),  This  program  sets  up  a  grid  of  the  kind  shown  in  Figure  4.1  for  a  cap  of 
size  CAPSIZ  ( in  degrees).  The  number  of  rings  in  the  cap  is  NP,  so,  with  the 
origin,  there  are  NP  +  1  weights  to  be  found.  The  listing  clearly  shows  the  various 
constants  used.  REGPAR  is  a  "regularization"  parameter,  chosen  here  very  small. 
This  number  is  added  to  the  main  diagonal  of  the  normal  matrix  C2Z  +  D  to  improve 
the  stability  of  the  solution.  NMOD  is  the  maximum  degree  and  order  in  the  refer¬ 
ence  model.  To  change  it  into  a  "perfect"  model,  a  statement  setting  all  values  in 
DVAR  to  zero  is  added  before  the  statement  "55  CONTINUE".  GNOISE  is  the  stan¬ 
dard  deviation  of  in  mgals.  The  solution  is  obtained  using  the  conjugate  gradi¬ 
ents  procedure  in  subroutine  CGRADS.  The  maximum  number  of  iterations  allowed, 
ITERMX,  is  the  number  of  unknown  (NP  +  1)  plus  ten.  However,  if  the  improvement 
in  the  rms  of  the  solution  from  iteration  to  iteration  is  less  than  one  part  in  a  million, 
the  procedure  terminates  then.  Program  RINGS  sets  up  the  normal  matrix  taking 
advantage,  as  far  as  possible,  of  the  various  symmetries  in  the  grid.  After  the 
solution  has  been  found,  it  multiplies  both  rms  and  weights  by  the  correction  factor 
1  +  kjE^  f4  (paragraph  (4.2)),  and  it  also  multiplies  (C2Z  +  D)_f  =  reconstructed 
right  hand  sides  of  normals,  to  compare  them  with  the  original  rms's,  showing 
whether  the  conjugate  gradient  procedure  has,  in  fact,  converged.  The  listing  of 
RINGS  is  complemented  by  those  of  the  subroutines  it  calls:  CGRADS,  MATVEC, 
LEGPOL,  COVAR,  and  function  F. 

A  A 

Finally,  subroutine  RINCOV,  used  to  determine  the  covariance  M{e  T(P4  )eT(P.)j 
between  prediction  errors,  based  on  (4.14)  and  the  repeated  use  of  expression  (4,17), 
is  listed  as  well.  The  main  array  and  variables  are  given  the  same  names  as  in 
RINGS  and  associated  subroutines,  with  the  exception  of  the  optimal  weights  vector, 
here  called  "F". 


SAMPLE  OUTPUT  OF  "RINGS" 
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COMPUTING  THE  COVARIANCE  BETWEEN  THE  ITH  AND  JTH  'RING  A' 
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